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FOREWORD 


The  work  described  in  this  report  was  performed  by  M.  Mahadeva 
Reddi  of  Conrad  Technologies ,  Inc.,  under  Contract  So.  S62269- 
86-C-0278  for  the  Naval  Air  Development  Center,  Warminster ,  PA. 
The  technical  monitor  was  Sr.  Lee  W.  Cause,  Advanced  Structures 
Technology,  NADC. 

Anticipating  that  the  readership  may  be  unfamiliar  with  some 
aspects  of  the  topics  covered,  an  attempt  has  been  made  to 
present  the  material  at  an  elementary  level.  While  this  may 
be  useful  to  some,  it  will  be  burdensome  to  others.  In  the 
latter  case,  skipping  the  section  may  provide  relief. 

The  author  is  particularly  grateful  to  Sr.  Cause  for  the 
valuable  technical  guidance  and  support  given  during  the  course 
of  this  work.  His  penchant  for  seeking  simplicity  when  faced 
by  seemingly  complex  situations  is  unparalleled. 

The  author  is  also  grateful  to  Dr.  Edmund  C.  Filetti  for  his 
assistance  in  performing  this  work. 
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1.  INTRODUCTION 

nte  general  procedures  for  structural  design,  analysis,  test 
and  data  requirenents  for  fixed  wing,  piloted  Navy  airplanes 
are  prescribed  in  the  llIL-A-8860  through  llIL-A-8870  series  of 
specifications  [1-12]*.  Within  the  spirit  and  intent  of  these 
specifications,  current  Navy  aircraft  design  practice  employs 
quasi-static  methodologies  to  determine  structural  response 
during  maneuvering  flight.  These  methodologies  are  based  on 
the  observation  that  during  a  maneuver  the  apparent  weight  of 
everything  in  the  aircraft  is  changed  by  a  load  factor,  N2, 
which  is  the  ratio  of  the  acceleration  of  the  aircraft  at  any 
instant  to  that  due  to  gravity.  The  load  factors  for  various 
flight  regimes  are  presented  in  the  form  of  what  is  called  the 
V-n  diagram,  shown  in  Fig.  1.1,  for  symmetric  flight  maneuvers. 
The  shape  of  this  diagram  is  derived  from  the  specifications  of 
the  procuring  activity  and  airworthiness  requirements  as 
prescribed  by  the  appropriate  authorities.  The  closed  region 
in  the  diagram  represents  a  maneuvering  envelope  for  the 
various  speed  and  load  combinations  which  may  be  encountered 
by  a  given  aircraft  during  its  design  life. 

The  specifications  enable  the  use  of  quasi-static  analysis  by 
prescribing  the  use  of  load  factors  to  account  for  the  tran¬ 
sient  loads  developed  during  maneuvers.  To  determine  the 
maneuver  load  factors,  an  initial  load  factor  for  each 
maneuver  is  prescribed  along  with  kinematic,  and  other 
parameters  of  the  maneuver  such  that  an  analysis  using  steady 
airloads  and  rigid  structure  assumptions  will  yield 
accelerations  at  the  center  of  gravity  of  the  aircraft  [2]. 

Numerals  in  square  brackets  designate  references 
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MIL-A.8861B(AS) 


MOTES: 

1.  JA  •  G8  •  vAlut  sptcifltd  In  columt  2  tnd  5,  tabli  I. 

2.  GC  •  valut  spactfltd  tn  colutnn  4,  tabit  I. 

3.  HO  >  KE  valut  sptcifltd  In  coluwis  3  and  6.  tabit  I. 

4.  OH  -  V  as  sptcifltd  In  mil-A.8860 

« 

5.  OG  -  V  as  sptcifltd  In  NIL-A-8860 

k 

6.  K  -  1.25  for  M-A  0.6 

-  1.0  for  M  >  1.0 

.  Cl. 625  -  (0.625  M)]  for  0.6<  N  <  1.0 

whtrt  M  is  tht  Mach  nu«btr  corrtsponding  to  tht  spttd  bting  considtrtd.  K  may 
bt  dtttrmintd  from  appllcabit  mind  tunnti  and  flight  ttst  data  acctptabit  to 
tht  procuring  activity.  This  dtttrmlnatlon  shall  Includt  considtratlon  of 
abruptness  of  tht  mantuvtr,  control  surfact  limitations,  Mach  numbtr,  thrust, 
canter  of  gravity  position,  external  stores  configuration,  maximum  safe  angle 
of  attack  as  limited  by  controllability,  limiting  buffet  loads,  and  other 
effects  which  can  be  shown  to  have  a  significant  bearing  on  the  maximum 
attainable  airplane  normal  force  coefficient  (C««  ) 

•I#  ■ 


Fig  1.1  V-n  Diagram  For  Symmetrical  Flight 
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Translated  into  load  factors,  these  accelerations  establish  the 
loads  to  be  applied  in  a  static  structural  analysis  to 
deternine  the  stresses  developed  during  the  naneuver.  By 
screening  the  results  of  the  aaneuver  transient  analysis, 
critical  "points-in-the-sky”  are  selected  for  detailed  stat'c 
structural  analysis  for  establishing  stresses. 

The  load  factor  approach  is  also  used  to  establish  the  fatigue 
design  service  loads  spectra,  nie  design  service  life  and 
design  usage  of  aircraft  is  specified  by  each  procuring 
activity  based  on  the  mission  requirements,  stipulated  in  terms 
of  total  flight  hours,  number  of  flights,  number  of  service 
years,  etc.  along  with  mission  profiles,  mission  nix  and  other 
special  requirements.  The  specifications  require  that  the 
design  service  loads  spectra  for  the  airframe  be  developed  for 
the  design  service  life  so  defined.  All  significant  sources  of 
repeated  loads  are  included  on  a  flight-by-flight  basis  to  form 
the  design  service  loads  spectra.  Aircraft  fatigue  analysis 
and  certification  testing  are  based  on  these  spectra,  the 
stresses  used  in  fatigue  analyses  again  being  derived  quasi- 
statically,  through  use  of  load  factors.  To  account  for  gust 
loads,  the  aircraft  is  presumed  to  encounter  gusts  in  both  the 
vertical  and  lateral  directions.  Hie  resulting  loads  are 
determined  by  the  discrete  gust  and  continuous  turbulence 
approaches  as  established  by  the  procuring  activity  [2]. 
Finally,  dynamic  analysis  is  required  for  assuring  freedom  from 
flutter,  divergence,  and  other  aeroelastic  instabilities  of  the 
aircraft  [9]. 

During  the  design  phases,  an  objective  in  the  iterative  loop 
of  design  modifications  is  to  achieve  an  economic  life  of  the 
aircraft  that  is  greater  than  the  design  service  life  when  the 
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aircraft  is  subjected  to  the  design  service  loads  spectra. 
Economic  life  of  the  aircraft  is  assumed  to  have  been  attained 
"with  the  occurence  of  widespread  damage  which  is  uneconomical 
to  repair  and,  if  not  repaired,  could  cause  functional 

problems  affecting  operational  readiness _ "[6].  To  achieve 

the  objective,  analytical  and  experimental  work  specified 
in  [6]  has  to  be  performed  to  demonstrate  compliance. 

Potential  Problem  Areas 

The  Navy  recognizes  that  the  use  of  accelerations  at  the  C.G. 
as  an  indicator  of  the  prevailii^  load  conditions,  is  valid 
only  for  structures  which  are  effectively  rigid.  While 
transport  aircraft,  which  are  flexible,  are  not  expected  to 
meet  this,  fighter  aircraft  are  considered  to  be  rigid  enough 
that  the  assumption  is  presumed  to  be  valid. 

The  current  design  methodology  deals  with  each  potential 
maneuver  separately,  by  isolating  it  from  the  flight  history 
preceding  it.  Although  an  initial  load  factor  is  prescribed 
for  maneuver  analysis,  whether  it  accounts  for  the  residual 
vibratory  effects  of  all  possible  immediately  prior  maneuvers 
is  not  evident.  In  fact,  the  inability  to  correlate 
instrumentation  data  with  analytical  results  [25]  would 
strongly  suggest  that  the  actual  peak  stresses  during  a 
maneuver  nay  be  significantly  different  from  the  results  of 
quasi-static  analysis,  because  of  the  residual  effects  of  the 
previous  maneuver.  If  this  is  so,  then,  there  may  be  a 
significant  shortfall  in  the  design  service  life. 

In  studying  the  development  of  maneuver  load  spectra  from  VGH 
(Velocity,  Load  factor  and  Altitude)  data,  Lincoln  [29]  states 
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that  "in  the  application  of  the  mission  analysis  required  by 
MIL-A-008666A  (USAF)  to  fighter  and  attack  aircraft,  a  problem 
arises  in  the  selection  of  the  point  in  the  sky  (velocity, 
altitude, and  weight)  for  the  load  factor  spectrum  for  the 
combat  segment  of  the  mission.  It  can  be  shown  that  in  many 
cases  important  differences  in  the  spectrum  can  be  obtained 
from  two  ’reasonable’  point  selections.” 

He  states  further  that  "the  problem  has  been  particularly 
severe  on  some  existing  aircraft  in  that  a  ten  percent  shift 
in  the  stress  spectrum  can  produce  a  factor  of  two  change  in 
life.  Therefore,  when  it  is  considered  that  essentially  all  of 
the  fatigue  damage  for  fighter  and  attack  aircraft  is  done  in 
the  combat  segment,  this  part  of  the  mission  deserves  special 
attention."  Needless  to  say,  maneuvering  loads  can  be  quite 
severe  in  a  combat  mission. 

Also,  to  maximize  the  economic  life  of  operational  aircraft, 
a  fleet  maintenance  program,  similar  to  (or,  same  as)  that 
specified  in  1IIL-STD-1530A  [13]  may  be  put  into  effect.  The 
basic  elements  of  such  a  plan  will  include  an  aircraft  tracking 
program,  a  load  spectra  survey,  and  a  fleet  structural 
maintenance  program.  A  key  ingredient  here  is  the  ability  to 
predict  potential  damage  in  critical  areas  of  each  airframe 
based  on  actual  usage  data  provided  by  instrumentation. 

However,  critical  locations  in  aircraft  are  frequently 
inaccessible,  and  it  may  not  be  practical  to  instrument  such 
locations,  and  indirect  means,  such  as  transfer  functions,  may 
have  to  be  used  to  predict  damage  at  those  locations  from 
observations  made  elsewhere  [14,15].  It  does  not  appear  to  be 
likely,  however,  that  such  functions  can  be  synthesized  without 
actual  dynamic  analysis.  If  inaccurate  transfer  functions  are 
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put  in  place,  damage  assessments  may  be  inadequate  and 
maintenance  actions  involving  repair  or  replacement  may  not 
occur  at  optimum  times,  nius,  there  can  be  a  shortfall  in 
terms  of  operational  costs  as  well. 

Moreover,  with  the  trend  toward  increased  optimization  in 
aircraft  design  through  the  use  of  computer  aided  analyses, 
aircraft  structures  are  becoming  more  and  more  "fully- 
stressed".  The  second  order  stresses  resulting  from  hitherto 
ignored  dynamic  stresses,  may  thus  become  significant  over  the 
life  of  the  structure. 

The  main  objective  of  the  effort  described  in  this  report  is, 
therefore,  to  explore  the  dynamics  of  the  internal  structural 
load  distributions  in  aircraft  during  maneuvering  flight  to 
establish  what  physical  phenomenon  are  significant  and  the 
level  of  analytical  modeling  necessary  to  make  conservative 
predictions. 

The  exploratory  study  consists  of  a  symmetrical  flying  wing 
simulated  by  simple  analytical  models  based  on  the  vortex 
lattice  method  for  quasi-steady  aerodynamics  and  time  domain 
solution  of  the  equations  of  motion  derived  by  the  finite 
element  method. 
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Structural  Oynanics  with  Independent  Loads 

The  current  state-of-the-art  for  analyzing  the  dynamics  of 
linear  structures,  subjected  to  loads  which  are  independent  of 
the  structural  deformations,  is  well  advanced,  dynamic 
response  can  be  readily  obtained  with  as  much  computational 
accuracy  as  is  desired.  However,  the  physical  properties  of  a 
structure  and  the  loading  conditions  are  generally  known  only 
approximately.  Hence,  the  structural  idealization  and  the 
solution  procedure  to  be  used  should  be  prudently  formulated  to 
provide  only  a  comparable  level  of  accuracy.  Practical 
problems  in  structural  dynamics  range  from  those  that  require 
only  models  of  extreme  simplicity  with  a  few  degrees  of  freedom 
and  requiring  only  one  or  two  modes  in  approximating  the 
dynamic  response  to  highly  sophisticated  finite  element  models 
with  several  thousand  degrees  of  freedom  in  which  hundreds  of 
modes  may  participate  significantly  in  the  response.  To  deal 
effectively  with  this  wide  range  of  analytical  requirements,  a 
variety  of  procedures  have  been  developed  which  have  proved  to 
be  efficient  in  practice. 

In  general  these  are  based  on  two  approaches,  both  of  which 
assume  that  a  transient  analysis  can  be  accomplished  most 
efficiently  by  establishing  first  the  structural  discretization 
necessary  for  static  stress  analysis  purposes,  and  then 
reducing  the  number  of  degrees  of  freedom  significantly  before 
performing  the  dynamic  analysis.  Hie  two  differ  in  the  manner 
in  which  the  number  of  degrees  of  freedom  is  reduced,  however. 
The  simpler  is  based  on  the  assumption  that  inertia  forces  are 
associated  only  with  certain  selected  degrees  of  freedom  of  the 
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original  discretization;  the  remaining  degrees  of  freedom  are 
not  explicitly  involved  in  the  dynamic  analysis  and  can  be 
ignored  in  the  dynamic  formulation.  In  the  second  approach, 
the  number  of  dynamic  degrees  of  freedom  is  limited  by  assuming 
that  the  displacements  can  be  combined  in  selected  patterns 
such  that  their  amplitudes  become  generalized  coordinates  for 
the  dynamic  analysis.  Several  specific  techniques  based  on 
these  can  be  found  in  [16]. 

Structural  Dynamics  with  Aerodynamic  Loads 


So  far  we  have  considered  structural  systems  subjected  to 
loads  which  remain  unchanged  while  the  structure  deforms.  The 
aerodynamic  loads  encountered  in  aircraft  structural  analysis, 
however,  are  for  the  most  part  functions  of  the  deformations 
and  their  rates  of  change  with  time.  Consequently,  the  overall 
problem  is  non-linear. 

To  make  the  problem  tractable,  a  particular  equilibrium 
configuration  (like  IG  level  flight)  is  selected  for  analysis, 
and  the  problem  is  linearized  by  assuming  small  perturbations 
about  the  equilibrium  configuration. 

The  traditional  procedure  for  obtaining  the  dynamic  response  of 
an  unrestrained  airplane  begins  with  modeling  of  the  continuous 
structure  by  one  consisting  of  a  finite  number  of  degrees  of 
freedom  and  the  aerodynamic  field  by  a  discrete  representation. 
This  leads  to  a  set  of  simultaneous  ordinary  differential 
equations  in  which  time  is  the  independent  variable  and 
structural  displacements  are  the  dependent  variables.  The  right 
hand  sides  of  the  equation  set  are  the  applied  forces  which 
usually  consist  of  disturbances  from  external  sources  such  as 
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gusts  and  those  arising  from  aircraft  notion.  By  an  eigenvalue 
analysis  of  the  homogeneous  set  of  equations  (obtained  by 
setting  right  hand  sides  to  zero)  the  natural  node  shapes  and 
frequencies  of  the  unrestrained  aircraft  structure  are 
obtained.  These  are  then  used  to  transform  the  original 
dynamical  equations  into  a  system  formulated  in  terms  of  the 
natural  modes,  which  then  become  the  generalized  coordinates 
for  further  analysis,  lliis  procedure  is  known  commonly  as 
solution  in  the  frequency  domain. 

For  gust  and  flutter  analyses,  the  equations  of  motion  have 
classically  been  solved  in  the  frequency  domain.  The  unsteady 
aerodynamic  forces  are  represented  by  functions  of  frequency. 

In  the  past  several  years,  however,  use  of  transform  techniques 
has  become  popular.  For  example,  with  Laplace  transforms,  the 
aerodynamic  loads  are  represented  by  complex  functions  of  the 
Laplace  variable.  The  frequency  dependent  tabular  represent¬ 
ation  of  the  loads  is  approximated  by  use  of  rational  poly¬ 
nomials  in  terms  of  the  Laplace  variable,  evaluated  by  a  least 
squares  fit  at  several  frequencies  of  interest.  This  tech¬ 
nique,  commonly  called  the  ”s-plane  approximation",  enables 
transformation  of  the  equation  of  motion  into  a  linear-time- 
invariant  (LTD  state-space.  Several  computational  techniques 
[32,33]  can  be  readily  applied  in  the  time  domain  of  the  LTI 
state-space;  hence  its  popularity. 

Traditional  aircraft  analytical  practice  also  distinguishes 
between  static  and  dynamic  loading  as  follows:  in  static 
loading,  the  inertia  forces  arising  from  vibratory  acceleration 
or  deceleration  of  the  masses  in  the  system  are  deemed 
negligible  compared  with  the  externally  applied  forces,  whereas 
in  dynamic  loading  they  are  considered  to  be  important.  For 
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static  loading,  thus,  the  steady  state  aerodynamic  forces  are 
equilibrated  with  inertial  forces  of  the  aircraft  which  is 
assumed  to  be  rigid.  For  the  dynamic  case,  the  inertial  and 
elastic  forces,  together,  are  set  into  equilibrium  with  the 
external  forces. 

Thus,  the  static  analysis  may  introduce  errors  in  two  ways: 
Structural  deformations  siay  produce  additional  aerodynamic 
forces  which  (1)  may  affect  the  overall  response  of  the 
aircraft  and  (2)  induce  vibrations  having  significant  effect  on 
internal  stresses.  Our  interest  in  this  survey  was  to  review 
methods  suitable  for  analyzing  response  to  rapidly  applied  or 
varying  loads.  However,  a  moderately  extensive  review  of 
available  literature  failed  to  disclose  any  particularly  well 
suited  advanced  methods.  A  brief  summary  of  some  articles 
which  are  representative  of  the  currrent  trends  follow: 

The  early  work  of  Bisplinghoff  et  al  [17],  distinguishes  three 
methods  for  computing  the  dynamic  loads,  namely,  mode 
displacement  method,  mode  acceleration  method,  and  summation  of 
forces  method.  For  a  sufficient  number  of  modes,  all  three 
methods  yield  analytically  identical  results.  However  for 
computational  economy  the  number  of  modes  required  for 
obtaining  given  accuracy  should  be  as  small  as  possible. 

Through  an  example,  Bishiplinghoff  et  al  have  demonstrated  that 
the  mode  acceleration  method  converges  faster  than  the 
mode  displacement  method.  These  methodologies  still  form  a 
part  of  much  contemporary  work. 

Some  of  the  analytical  tools  used  at  Lockheed  to  design  trans¬ 
port  aircraft  to  survive  the  transient  environments  experi¬ 
enced  during  their  service  life  are  described  by  Wignot  [181. 
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Flight  transients  include  gusts,  steady  aaneuvers  and  abrupt 
maneuvers,  and  ground  transients  (which  include  landing, 
taxiing,  braking  and  turning).  All  such  analytical  methods  are 
based  on  the  rationale  that  the  time  rate  of  change  of  notion 
during  a  maneuver  is  slow  enough  that  the  maneuver  can  be 
treated  as  a  steady-state  aeroelastic  problem.  Brief 
descriptions  of  four  computer  programs,  for  steady  symmetric 
flight  pitch  maneuver,  transient  symmetric  flight  pitch 
maneuver,  lateral  maneuver,  and  lateral-  directional  maneuver 
are  given.  Tine  history  analyses  are  also  described  but  the 
computational  techniques  used  in  the  analyses  are  not 
elucidated.  The  paper  also  describes  in  a  general  way  program 
capabilities  for  analyzing  impulsive  loadings  such  as  hail  or 
bird  strike,  rapid  decompression,  and  abnormal  landing  gear 
configuration  landings. 

An  analytical  technique  for  determining  the  dynamic  response  of 
flexible  aircraft  is  described  by  Vogel[19].  The  excitations 
considered  are  gusts  and  arbitrary  maneuvers.  The  solution  is 
based  on  an  "analytic  representation  of  the  admittance 
functions  by  partial  fractions".  The  time  response  frequency 
dependence  of  the  aerodynamic  forces  induced  by  motion  are 
approximated  to  allow  analytical  solutions.  "The  results 
are  shown  to  be  as  accurate,  but  more  economical,  than  more 
elaborate  numerical  Fourier  methods.  Elasticity  is  introduced 
through  eigenf requeue ies  and  generalized  masses.  Several 
examples  are  given  shich  also  include  the  influence  of  control 
systems  and  of  elasticity  on  loads  and  accelerations".  The 
technique  is  claimed  to  be  an  inexpensive  approximating  method 
for  all  dynamic  problems  in  which  "elasticity,  control  systems 
and  excitation  time  lag  effects  play  a  role". 
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Msirtin  [20]  describes  a  aatrix  load  analysis  nethod  for 
flexible  aircraft  structures.  The  overall  objective  is  to 
account  for  the  effect  of  flexibility  on  load  distribution. 

Ihe  aerodsmaaic  loads  are  conputed  on  the  basis  that  the  final 
load  distribution  equals  the  initial  distribution  for  a  rigid 
aircraft  plus  the  change  in  distribution  due  to  its  structural 
flexibility.  Dsmanics  are  not  included  in  the  analysis. 

Pototzky  and  Perry  [21]  review  existing  techniques  for 
calculating  dynamic  loads  for  flexible  airplanes  and  present  a 
new  technique  based  on  the  suamation-of-forces  method  of 
writing  dynamic  loads  equations.  This  paper  is  noteworthy  for 
its  concise  review  of  the  state-of-the-art  in  dynamic  loads 
analyses.  They  also  demonstrate  that  the  summation  of  forces 
method  converges  faster  than  the  mode  displacement  method. 

A  computer  program  for  flexible  aircraft  flight  dynamic  loads 
analyses  with  active  controls,  entitled  DYIOFIEX,  is  described 
by  Perry  III  et  al[22].  The  dynamical  equations  are  formulated 
with  the  assumptions  that  all  notions  are  small  and  that  the 
aircraft  is  in  trimmed,  unaccelerated,  straight  and  level 
flight.  Time  histories  of  dynamic  loads  due  to  discrete 
excitations  (gust  and  control  surface)  are  computed  by  Fast 
Fourier  Transform  techniques  using  the  Cooley-Tookey 
algorithm  [23].  Both  steady  and  unsteady  aerodjmanics  are 
included.  Structural  data  in  the  form  of  mode  shapes, 
generalized  mass  and  stiffness  matrices  generated  by  external 
programs,  finite  element  or  otherwise,  are  required  as  input 
to  the  program. 


Vakhitov  et  al  [24]  developed  a  numerical  procedure  for 
determining  stresses  in  a  wing  under  transient  loading.  The 
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aerodynamic  load  is  determined  by  a  discrete  vortex  method. 

The  djmamical  equations  for  the  flexible  structure  are 
integrated  forward  in  time  using  a  finite  difference  approach. 
The  response  of  a  wing  subjected  to  gust  loading  is  given  as 
an  example.  This  article  is  analytically  the  closest  to  the 
approach  used  in  the  current  analysis  as  described  in  this 
report. 
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3.  SYMMETRIC  WING  MODEL 


We  know  fron  ordinary  experience  that  if  a  load  is  applied  to 
a  flexible  structure  slowly  enough,  then  the  structure  will  not 
exhibit  an  overswing.  To  S3mthesize  a  uodel  for  the  purposes 
of  the  present  study,  it  is  instructive  to  examine  how  slowly 
the  load  has  to  be  applied  in  order  that  the  overswing  will  be 
small.  For  this  purpose  we  consider  an  undamped,  single  degree 
of  freedom  oscillator  subjected  to  loading  which  increases  to  a 
maximum  value  of  F  in  tine  r  and  remains  constant  thereafter. 
During  the  interval  0  <  t  <  r,  the  dynamical  equation  for  the 
single  degree  of  freedom  system  is 

m^+kx  =  F|-  3.1 

dt’  ^ 

the  general  solution  of  which  is 


X  a  ^  +  A  sin  wt  +  B  cos 


where, 


Assuming  that  the  system  is  at  rest  at  time  t  s  0, 


the  solution  is 


The  elastic  force,  -kx,  can  now  be  expressed  as  a  function  of 
tine,  as  follows 


NAOC  88014-60 

Symmetric  Wing  Model 


-to  .  -  f  t  ♦  f  3.6 

where  the  firet  tera  on  the  right  hand  tide  balances  the 
applied  force  and  the  second  tera  balances  the  inertia  force. 
As  tiae  increases  froa  zero  to  r,  the  applied  force  also 
increases  but  the  inertia  force,  however,  oscillates  between 
±  Therefore,  if  u  is  large  enough,  then  the  inertia 

force  will  becoae  negligible  in  coaparison  with  the  aaziaua 
value  of  the  applied  force  F.  Thus,  as  the  load  application 
tiae  becoaes  greater  in  coaparison  with  the  period  of  natural 
vibration,  then  the  overswing  due  to  inertia  force  becoaes 
smaller. 


For  t  >  r,  the  solution  is  given  hy 

*  »  f  II  ♦  j3^-sin  ui  *  sin  ((/(t  -  r)}]  3.6 

and  the  elastic  force  is 

-kx  =  -F  tl  ♦  j3^-sin  wt  +  sin  w(t  -  r)}]  3.7 

The  fractional  overswing  of  the  elastic  force  over  the  applied 
force  is  thus  of  the  order  of  l/urr  for  the  undaaped,  single 
degree  of  freedoa  oscillator.  For  systeas  with  aultiple 
degrees  of  freedoa,  by  extension,  we  aay  conclude  that  the  rate 
of  load  application,  together  with  the  lower  natural 
frequencies,  will  affect  the  aaount  of  overswing  siailarly. 


Our  aain  guidelines  in  synthesizing  a  aodel  for  investigating 
the  severity  of  aaneuver  transients  are  siaplicity  and 
conservatisa.  By  presuaing  that  siaple  aaneuvers  will  need 
only  siaple  aodels  to  siaulate  thea,  we  define  a  aaneuver  for 
synthesizing  the  aodel  as  follows:  An  aircraft  flying  with  an 
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initial  load  factor  enters  soae  Maneuver  in  a  tine  r  such  that 
after  all  transients  subside,  a  different  load  factor  is 
attained.  Our  objective  is  to  detemine  what  the  naxinun 
overswing  is  by  the  tine  steady  state  is  reached.  We  also 
presune  that  the  Maneuver  is  synetric  so  that  the  size  of  the 
coaputational  burden  nay  be  cut  in  half  or  aore. 

The  sinplest  nodel  which  can  be  subjected  to  such  a  naneuver 
is  a  subsonic  flying  wing.  Its  structure  can  be  nodeled  as  a 
bean  and  its  aerodynanics  can  be  detemined  by  vortex  lattice 
Method.  Fig.  3.1  shows  typical  layouts  of  the  structural 
nodes  and  the  vortex  lattice  panels. 


Vortex  Lattice  Panels 


Finite  Clenent  Nodes 


Fig  3.1  Typical  Finite  Eleaent  and  VUI  Panel  Layout 
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4.  VORTEX  LATTICE  METHOD 


After  a  brief  review  of  subsonic  airfoil  theory,  the  numerical 
approach  for  arbitrary,  finite  planar  wings  using  vortex  panels 
is  described.  The  development  follows  that  of  Bertin  and  Smith 
[31]  closely.  (Aerodjmamicists  may  safely  skip  this  section.) 


Elementary  aerodynamics  teaches  us  that  superposition  of  three 
elementary  flows  -  uniform  flow,  source  flow  and  doublet  flow  - 
can  be  used  to  obtain  nonlifting  flow  around  several  body 
shapes  such  as  cylinders,  ellipsoids  etc.  and  that  a  fourth 
elementary  flow  -  vortex  flow  -  is  necessary  to  obtain  lifting 
flow.  In  vortex  flow,  all  the  streamlines  are  concentric 
circles  with  the  tangential  velocity,  V^,  varying  inversely 
with  radius,  the  radial  velocity  being  zero.  Thus  we  have, 


=  0  4.1 

4.2 

The  constant  in  Eq.  4.2  nay  be  evaluated  from  the  definition  of 
circulation,  F,  around  a  given  streamline  of  radius  r. 


r  =  d)  V  .  ds  =  -V/ 


2wr 


4.3 


or 


Va  = 


4.4 


2Tr 

The  minus  sign  arises  from  the  definition  that  positive  circu¬ 
lation  is  clockwise,  and  consequently,  a  vortex  of  positive 
strength  rotates  clockwise. 


The  Kutta-Joukowski  theorem  states  that  the  lift  per  unit  span, 
1 ,  on  a  two  dimensional  cross-section  of  a  body  is  directly 
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proportional  to  the  circulation  around  the  body.  Thus, 

1  -  P«V«r  4.5 

where.  is  the  uniform  flow  velocity  at  infinity,  is  the 
fluid  density  and  T  is  the  circulation  around  the  body.  Ihe 
circulation  may  be  calculated  over  any  closed  path  that 
encloses  the  body. 

Tht  idea  of  a  point  vortex  in  two  dimensional  space  can  be 
extended  to  include  a  line  vortex,  and  indeed  even  a  sheet 
vortex,  when  three  dimensional  space  is  considered;  the  vortex 
strength  7  is  then  simply  a  spatial  function.  To  compute  lift 
on  an  arbitrary  body,  a  vortex  sheet  is  placed  on  the  body’s 
surface,  the  vortex  strength  being  such  that  the  induced 
velocity  field,  when  combined  with  the  uniform  flow  makes  the 
vortex  sheet  (consequently  the  body  surface)  a  streamline  of 
the  flow.  The  circulation  is  in  turn  found  from 


where  ds  is  the  appropriate  length  or  area  element,  the 
integral  being  taken  around  the  complete  surface  of  the  body. 
The  lift  is  found  as  usual  from  eq.  4.5.  However,  for 
arbitrary  body  geometries,  a  general  analytical  solution  for 
7  s  7(s)  has  not  been  found.  Rather,  it  has  to  be  found  by 
numerical  methods. 

The  considerations  discussed  so  far  do  not  by  themselves  assure 
uniqueness  of  the  solution  for  7(0).  Depending  on  the  circu¬ 
lation,  r,  an  infinite  number  of  potential  flow  solutions  are 
possible.  To  select  the  proper  one.  an  additional  condition 


NADC  88014-60 

Vortex  Lattice  Method 

has  to  be  inposed.  This  is  given  by  the  Kutta  condition,  which 
states  that  the  value  of  F  is  such  that  the  flow  leaves  the 
trailing  edge  smoothly.  More  precisely,  if  the  trailing  edge 
has  a  finite  angle,  then  the  trailing  edge  is  a  stagnation 
point;  if  the  edge  is  cusped,  then  the  flow  velocities  leaving 
the  top  and  bottom  surfaces  are  equal  at  the  edge  in  magnitude 
and  direction. 

When  finite  wings  are  considered,  because  of  the  higher 
pressure  at  the  bottom,  a  flow  from  the  bottom  to  the  top 
develops  at  the  wing  tips.  Also,  spanwise  pressure  gradients 
are  such  that  secondary  flows  develop  from  the  tip  to  the  root 
on  the  top  and  from  the  root  to  the  tip  on  the  bottom  surfaces. 
These  secondary  flows,  in  particular  that  around  the  wing  tips, 
create  vortices  around  the  tips  which  trail  downstream  of  the 
wing  inducing  a  small  velocity  component  in  the  downward 
direction  called  a  downwash.  The  presence  of  downwash  modifies 
the  relative  wind  direction  seen  by  the  wing  such  that  the 
angle  of  attack  and  the  lift  direction  are  modified.  To  compute 
the  downwash,  we  use  the  Biot-Savart  Law  and  the  Helmholtz 
theorem,  which  state  that  the  strength  of  a  vortex  filament  is 
constant  along  its  length  and  that  a  vortex  filament  cannot  end 
in  the  fluid. 

Numerous  computational  schemes  have  been  developed  to  predict 
the  flow  around  a  thin  wing  flying  with  a  small  angle  of  attack 
so  that  steady-state,  inviscid,  irrotational  and  incompressible 
assumptions  may  be  invoked.  In  the  vortex  lattice  method  (VLM), 
the  continuous  distribution  of  a  vortex  sheet  is  approximated 
by  emplacing  a  number  of  discrete  horseshoe  vortices  in  trape¬ 
zoidal  panels  as  shown  in  Fig.  3.1.  The  bound  vortex  in  each 
panel  is  placed  in  such  a  way  that  it  is  coincident  with  the 
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quarter  chord  line  which  is  aligned  with  the  local  sweepback 
angle.  In  a  rigorous  analysis  the  vortex  lattice  panels  would 
be  located  on  the  mean  camber  surface  of  the  wing  and  the 
trailing  vortices  would  follow  a  curved  path  as  they  leave. 

For  our  present  purposes,  since  a  thin  flat  wing  is  judged  to 
be  adequate,  suitable  accuracy  can  be  obtained  by  using  linear 
theory  in  which  straight-line  trailing  vortices  aligned 
parallel  to  the  vehicle  axis  extend  downstream  to  infinity. 
This  orientation  of  the  trailing  vortices  yields  a  simpler 
computational  procedure  than  one  aligned  parallel  to  the  free 
stream  while  giving  similar  accuracy.  Also,  various  geometric 
coefficients  remain  independent  of  the  angle  of  attack.  For 
our  transient  dynamic  analysis,  this  is  a  desirable  feature 
because  it  speeds  up  the  computation  by  allowing  decomposition 
of  a  matrix  only  once.  Subsequent  continued  solutions  then 
require  only  changing  right  hand  sides. 

The  flow  field  induced  by  a  horseshoe  vortex  is  obtained  by 
superposition  of  the  flows  due  to  each  of  the  three  straight 
line  segments  determined  separately.  The  Biot-Savart  law 
enables  us  to  compute  the  flow  dV  induced  by  a  vortex  filament 
of  strength  F  and  length  dL  as 

dV  _  r  tdi  X  R)  4.7 

45rr^ 

In  a  typical  wing  panel,  say  the  nth,  let  the  beginning  and 
end  points  of  the  bound  vortex  be  given  by  yjn* 

(Xjni  Yjn*  *jn^’  integrated  over  the  bound 

vortex  and  the  two  trailing  vortices  extending  to  infinity  to 
obtain  the  velocity  induced  at  some  point  (x,y,z}. 

The  total  velocity  Vg|  induced  at  the  mth  control  point  by 
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the  nth  panel  may  be  stated  as 
'^m.n  “  *^,n^n 

where  the  influence  coefficient  depends  on  the  geometry 
of  the  nth  vortex  and  its  distance  from  the  mth  panel.  The 
velocities  induced  by  all  the  vortices  may  be  superposed  to 
give  the  velocity  induced  at  the  mth  control  point: 

2N 

n=l 

There  are  2N  such  equations,  one  for  each  panel.  If  the 
strengths  of  the  2N  vortices  are  known,  then  it  is  possible  to 
compute  the  induced  velocity  at  any  point  in  space.  To 
determine  the  strengths  we  use  the  boundary  condition  that  the 
flow  is  tangent  to  the  panel  at  the  control  point. 

For  a  planar  wing  in  the  xy  plane,  with  the  trailing  vortices 
parallel  to  the  x-axis,  all  three  components  of  the  vortex  of 
the  mth  panel  induce  a  velocity  at  the  control  point  of  the 
nth  panel,  parallel  to  the  z-axis  (  a  downwash).  Thus,  eq. 

4.9  can  be  written  as 


2N 

■  nil 


4.10 


where  w„  is  the  downwash  at  the  control  point  of  the  mth  panel 
m 

and  w_  -  is  the  downwash  induced  at  the  control  point  of  the 

iD,n 

mth  panel  by  the  vortices  of  the  nth  panel,  (derived  in  Bertin 
and  Smith  [31]),  as  follows: 
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m.n 


47r 


m  in 


A  / (*  ~x  )*  +  (y  -y  )' 
^  m  2n  m  2n 


4.11 


y  _  y 

■'in  •’m 


1  + 


X  -X. 

m  in 


in  in  m  in 


^2n 


1  + 


* 

m  2n 


m  2n 


m  •^2n 


The  component  of  the  combined  flow  normal  to  the  wing  is 
required  to  be  zero.  The  free-stream  velocity  component  normal 
to  the  wing  is  V^sin  a,  where  a  is  the  angle  of  attack.  The 
oscillatory  wing  velocity  normal  to  the  wing  is  W^^j.  Combining 
them , 


Wm  sin  a  =  0 

m  In  QD 


4.12 


Applying  Eq.  4.12,  the  strength  of  the  vortices  may  be 
determined  and  the  lift  can  then  be  calculated  from  Eq.  4.5. 


Unsteady  Aerodynamics 


A  parameter  that  characterizes  the  effects  of  unsteady  aero¬ 
dynamic  effects  is  the  reduced  frequency,  k,  given  by 
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where  c  is  the  semi-chord  and  is  the  radian  frequency. 

The  parameter  k  is  a  measure  of  the  importance  of  unsteady 
effects  in  the  flow,  and  is  2w  times  the  ratio  of  the  time  it 
takes  the  wing  to  travel  one  half  chord  length  to  the  time  for 
one  complete  oscillation.  If  k  is  small,  the  oscillatory 
effects  of  the  wing  motion  are  small  compared  to  the  forward 
motion  and  the  fluid  dynamics  can  be  represented  by  quasi¬ 
steady  modeling.  Thus,  at  any  given  instant,  the  wing  may  be 
assumed  to  be  fixed  and  steady  state  aerodynamics  may  be  used 
to  calculate  fluid  forces.  The  wing  position  is  then  advanced 
and  the  computation  is  repeated  for  the  next  instant. 

The  value  of  k  for  which  unsteady  effects  become  important 
may  be  judged  from  the  observation  that  for  k=0.2,  with  two- 
dimensional  or  high  aspect  ratio  wings,  unsteady  effects 
become  fully  evident.  For  the  model  being  considered  here, 
k  ranges  from  .12  at  the  tip  where  the  vibratory  amplitude  is 
a  maximum  to  .47  at  the  root  where  the  amplitude  is  a  minimum, 
that  at  the  mean  chord  being  .295.  Thus,  the  effect  of 
unsteady  aerodynamics  can  be  fully  developed  in  the  model,  but 
because  of  its  computational  speed,  quasi-steady  modeling  was 
selected  as  the  preferred  alternative  to  incorporating  the 
vortex  lattice  method  for  general,  unsteady  aerodynamics  as 
described  by  Konstadinopoulos  et  al  [34]. 
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5.  FINITE  ELEMENT  FORMULATION 


Finite  element  methodology  for  structural  applications  is 
reviewed  briefly.  (Structural  dynamicists  may  safely  skip  this 
section. ) 

In  problem  solving,  if  the  original  problem  is  too  difficult 
for  a  direct  solution,  then  a  strategy  that  often  works  is  to 
reduce  it  to  a  set  of  smaller  problems  which  may  be  easier  to 
solve.  One  may  view  the  application  of  the  finite  element 
method  to  structural  problems  along  these  lines  by  visualizing 
a  structure  as  an  assemblage  of  smaller  structures,  inter¬ 
connected  at  points  called  nodes.  Obviously  the  smaller 
structures  must  be  such  that  the  problems  they  present  are 
tractable.  Thus,  the  wing  structure  of  an  airplane  may  be 
viewed,  for  the  present  purposes,  as  an  assemblage  of  bean 
elements,  with  the  proviso  that  for  each  element,  we  know  how 
to  predict  the  structural  behaviour.  More  specifically,  if 
forces  and  moments  are  applied  at  the  nodes  of  each  element, 
then  we  know  what  the  resulting  deflections  are. 

The  formal  and  elegant  way  of  calculating  the  element 
properties  is  through  variational  principles.  However,  for  the 
present  purposes,  it  will  suffice  for  us  to  review  rather 
briefly  the  elementary  direct  method  which  results  naturally 
from  the  basic  techniques  of  structural  analysis.  In  cases 
where  the  structure  may  be  idealized  by  an  assemblage  of  one 
dimensional  elements  alone,  the  direct  method  is  just  as 
powerful  as  the  variational  method. 

The  basis  of  the  direct  method  is  to  find  stiffness  influence 
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coefficients  such  that  th«  generalized  force  at  one  point  of 
the  structure,  arising  fron  a  unit  generalized  displacenent  at 
the  same,  or  a  different  point  in  the  structure,  can  be 
•  determined.  For  a  bean  element,  for  example,  if  all  the  nodal 

displacements  (including  rotations)  except  for  one  are  set  to 
zero  and  the  free  one  set  to  unity,  then  the  nodal  forces  (and 
moments)  necessary  to  maintain  such  deformation  of  the  element 
constitute  a  colusm  of  what  is  termed  the  stiffness  matrix, 
the  column  number  corresponding  to  the  generalized  displacement 
that  is  set  to  unity.  To  illustrate,  if  a  bean  element  has 
four  generalized  degrees  of  freedom  and  the  second  is  set  to 
unity  while  the  others  are  kept  at  zero,  the  required 
generalized  forces  constitute  the  second  column  of  the  matrix. 
Thus,  in  the  following  static  equilibrium  equation 

Kq  =  f  5.1 

where  R  is  the  stiffness  matrix  and  q  and  f  are  vectors  of 
generalized  displacements  and  forces  respectively,  the 
elements  of  K,  q  and  f  nay  be  explicitly  written  out  as 


Tr  Ir  V  Ir  ' 

•‘ll  *^12  *13  *14 
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k  k  k  k 

Rjj  Rjj  *23 

<»2 

^2 

k  k  k  k 
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f3 

k  k  k  k 

L  41  42  *43  *44j 
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5.2 


For  beam  elements,  the  generalized  forces  are  the  force  and 
moment  applied  at  each  node  and  the  corresponding  generalized 
displacements  are  the  lateral  displacements  and  rotations 
respectively. 

Under  the  assumptions  of  plane  stress,  constant  flexural 
rigidity  (El)  and  no  shear  deformation,  by  considering  four 
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diiplaceaent  vectors  in  which  all  displaceaents  except  for  one 
are  set  to  zero,  the  coefficients  in  Eq.  5.2  say  be  evaluated 
froB  elementary  structural  nechanics.  Eq.  5.3  shows  the 
derivation  for  two  coluauis.  The  eleaents  of  the  stiffness 
Batrix  are  then  as  follows: 


where. 
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5.3 


5.4 


The  remaining  elements  follow  from  sysusetry,  namely, 


5.5 


System  Assembly 


To  represent  the  structure,  the  individual  finite  elements 
have  to  be  assembled  into  the  overall  system  stiffness  matrix. 
Referring  to  Fig  5.1,  two  adjoining  beam  elements  are  shown 
prior  to  assembly  and  after  assembly.  At  the  node  joining  the 
two  elements,  the  applied  forces  are  the  sum  of  those  applied 
to  the  individual  elements.  Also,  the  displacements  at  the 
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no4t  art  tha  taae  as  thosa  for  tha  individual  alanantt,  because 
of  physical  continuity  raquiraaants.  Consaquantly,  tha 
ovarall  static  aquilbriua  aquation  for  tha  two  assaablad 
alaaants  can  ba  fomulatad  as  follows: 


K 

1 

■>; 

ft 

«n**n*i 

«!*■ 

'  m  • 

fj  ♦ 

fj  4.  f5^» 

*n+i 

_n+i 

1  **4  J 

fn<fi 

fn-fi 

L  4 

Fig. 5.1  Baaa  Elaaants 


Tlia  systaa  stiffnass  aatrix  is  populatad  with  tarns  fron  each 
of  tha  alanant  stiffnass  natricas  as  shown  in  Eq.  5.6.  The 
partitioned  areas  labeled  and  contain  tarns  froa  finite 
elaaants  n  and  n+1,  respectively;  their  intersection,  labeled 
contains  tarns  which  are  tha  suas  of  tha  respective 
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terms  from  the  individual  stiffness  matrices,  llius,  the  system 
stiffness  matrix  may  be  assembled  by  repeated  application  of 
the  foregoing  algorithm  to  all  the  elements  in  the  structural 
assemblage. 

Mass  Matrix 


Prediction  of  the  dynamic  behaviour  of  a  structure  requires  not 
only  its  stiffness  properties  but  also  the  inertial  properties. 
For  finite  element  analysis,  the  inertial  properties  are  cast 
into  the  form  of  a  mass  matrix,  the  two  basic  types  being  the 
consistent  mass  matrix  and  the  lumped  mass  matrix.  Whereas  the 
former  is  based  on  variational  principles,  the  latter  is 
obtained  more  directly  by  lumping  the  mass  of  an  element  at  its 
nodes  by  assuming  that  the  mass  adjacent  to  a  node  is 
concentrated  at  the  node.  Ihus,  for  a  uniform  bar  element,  the 
isass  of  the  bar  is  equally  apportioned  between  the  two  nodes. 
For  our  present  purposes,  the  lumped  mass  approach  is  judged  to 
be  more  than  adequate. 

Reduction  of  Unloaded  Degrees  of  Freedom 


In  modeling  an  airplane  wing  with  beam  elements,  if  the  effects 
of  rotatory  inertia  are  deemed  to  be  negligible,  as  is  the  case 
in  the  present  application,  and  no  external  moments  are  applied 
to  the  nodes,  then  in  the  equilibrium  equations,  the  right  hand 
sides  for  the  corresponding  degrees  of  freedom  are  zero.  This 
fact  may  be  taken  advantage  of  to  lessen  the  computational  bur¬ 
den  by  reducing  the  total  number  of  degrees  of  freedom  involved 
in  the  dynamic  response  computation. 

Thus,  if  the  equilibrium  equations  are  rearranged  and 
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partitioned  as  follows. 


where  now  represents  the  vector  of  rotational  degrees  of 
freedom  which  are  not  subject  to  any  externally  applied 
moments,  one  can  solve  for 

^2  “  “  *22  *21^1 

Substituting  into  Eq.  5.7  for  Q,. 

[Kjj-  KjjKjjKjj]  Qj  =  Fj  5.9 

where  the  quantity  in  the  square  brackets  is  the  effective 
stiffness  matrix  for  the  translational  degrees  of  freedom. 

This  process  is  called  static  condensation  or  reduction. 

The  reduction  procedure  as  described,  while  formally  correct, 
is  nevertheless  computationally  inefficient.  The  more 
efficient  way  is  to  carry  out  a  symmetric  Gaussian  elimination 
backwards.  A  Fortran  version  of  this  procedure  may  be  found  in 
the  Appendix. 
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6.  DYNAMICAL  EQUATIONS 


A  Mneuvering  aircraft  in  flight  ia  a  three  dimensional, 
flexible  body,  free  to  translate  and  rotate  in  space.  It  is 
subject  to  aerodynamic  forces  sdiich  are  themselves  functions 
of  the  deformations  of  the  aircraft  structure  because  of  its 
flexibility. 

In  a  Newtonian  reference  frame,  the  dynamical  equations  for 
our  flying  wing  are  given  by 

M^  +  C^  +  KW  =  F  6.1 

dt*  dt 

where,  U  is  a  displacemnt  vector.  N,C.K  are  the  mass,  damping 
and  stiffness  matrices,  respectively,  and  F  is  the  applied 
force  vector.  In  general,  the  force  vector  F  can  be  a 
non-linear  function  of  the  displacement  vector  N  and  as  such 
the  system  of  equations  represented  by  Eq.  6.1  can  be 
non-linear.  Also,  the  initial  conditions  can  be  arbitrary. 

The  direct  method  for  solving  Eq.  6.1  is  to  integrate  it  in 
time.  However,  for  systems  involving  large  number  of  degrees 
of  freedom,  this  can  be  expensive.  It  is  possible  to  reduce 
the  computational  burden  if  the  number  of  modes  can  be 
reduced.  This  is  possible  because  generally  only  the  lower 
modes  participate  significantly  in  the  motion. 

Structures  which  are  unconstrained  against  rigid  body  motions, 
such  as  aircraft  or  missiles,  require  special  analytical 
treatment  because  the  stiffness  matrices  are  singular  with  the 
eigenvalues  corresponding  to  the  rigid  body  modes  being  zero. 
In  analyses  formulated  with  intermediate  results  obtained  in 
the  frequency  domain,  the  stiffness  matrix  is  swept  to  remove 
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the  offending  nodes.  Similar ly,  with  narching  schenes,  such 
as  Runge-Rutta,  the  rigid  body  notions  soon  become  so  much 
larger  than  the  defornations  in  the  structure  that  computa¬ 
tional  accuracy  of  the  defornations  is  totally  lost  and  the 
solution  becomes  unstable.  One  of  the  ways  in  which  this  can 
be  avoided  is  to  employ  transformation  of  coordinates  to  an 
accelerated  frame  by  eliminating  rigid  body  notions  at  each 
step  of  the  marching  scheme. 

For  this  purpose,  we  will  employ  a  transformation  such  that 
the  elastic  deformations  will  be  determined  in  a  coordinate 
system  fixed  at  the  C.G.  Thus,  in  an  inertial  reference 
frame,  let  R  be  the  position  vector  of  the  C.G.  Then, 

U  =  R  +  U  6.2 

where,  I)  is  the  vector  of  elastic  defornations  with  reference 
to  the  C.G.  Substituting  Eq.  6.2  into  6.1,  we  have 

M  ^  +  C  ^  ♦  K  U  =  F’-  M  6.3 

where  the  elastic  and  internal  damping  forces  for  a  rigid 
body  displacement  have  been  set  to  zero.  Also,  F’  now 
includes  aerodynamic  contributions  arising  from  wing 
vibrations . 

The  initial  conditions  corresponding  to  steady  level  flight 
under  IG  loading  will  ensure  a  minimum  starting  transient. 

For  this  purpose  we  first  obtain  a  static  solution  under 
Ig  loading.  These  initial  conditions  nay  be  stated  as: 
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7.  THE  RUWGE-KUTTA  METWOD  FOR  DIFFERENTIAL  EQUATIWS 
OF  THE  SEOCWD  ORDER 


nie  Runge-Kutta  method  is  a  popular  technique  for  integrating 
ordinary  differential  equations.  The  general  formulae  for  the 
algorithm  specialized  to  various  orders  may  be  found  in  [261. 

We  detail  here  the  special  case  of  second  order  equations  which 
is  particularly  relevant  to  structural  problems.  This  approach 
is  favored  because  of  the  symmetry  and  simplicity  of  the 
coefficients  and  an  accuracy  of  the  fourth  order. 

Given  the  vector  of  dependent  variables  Y  and  the  independent 
variable  x,  let  the  differential  equation  be  of  the  form 

y”  =  f(x.  Y.  y’)  7.1 

where  the  prime  denotes  differentiation  with  respect  to  the 
independent  variable  x.  The  integration  is  subject  to  the 
initial  conditions 


Yo  =  Y  (Xj,) 


For  computational  ease,  it  is  convenient  to  work  with  the 
quantities 


V  =  h  y’  7.3 

kj,  =  f(x,  Y,  V/h) 


where,  h  is  the  increment  of  the  independent  variable  x  for 
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the  integration  step,  and  =  1,  2,  3,  or  4. 

The  appropriate  formulae  for  the  intermediate  computations  are 
given  in  Table  7.1. 

Table  7.1  Runge-Kutta  Formulae  for  y’’=  f(x,  Y,  y’ ) 


i 


If  the  function  f  does  not  contain  Y  explicitly,  then  kj  s 
and  consequently  there  will  be  one  less  row  to  be  computed. 

Collatz  [27]  offers  a  rule  of  thumb  for  estimating  the  size  of 
the  increment  h:  ”  ...  h  should  be  such  that  kj  and  kg  are 
approximately  equal;  to  be  more  specific,  the  difference 
between  kg  and  kg  should  not  exceed  in  magnitude  a  few  percent 
of  the  difference  between  k^and  kg,  otherwise  a  smaller  step 
should  be  taken.  The  ratio  (kg-  kg)/(kj-  kg)  is  a  measure  of 
sensitiveness. . .” 
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8.  CXMPUTER  PROGRAM  DESCRIPTI(»4 


Hie  main  program  for  performing  the  maneuver  dynamic  analysis 
is  called  MANDYN.  Written  in  Fortran  77,  it  was  compiled  and 
executed  with  Microsoft  Fortran,  Vers.  4.0  on  an  IBM  PC-AT, 
with  a  co-processor,  and  a  20M  hard  disk.  Typical  run  times 
were  about  5  minutes  for  2000  time  steps  for  a  model  comprising 
9  each  of  beam  finite  elements  and  VLM  panels.  PC-DOS,  Vers. 
3.2  was  the  operating  system.  A  source  listing  of  the  program 
can  be  found  in  the  Appendix.  Input  to  the  program  is  menu 
driven. 

The  main  module  MANDYN  provides  user  interface,  initializes 
the  program,  marches  through  the  time  steps  and  creates  files 
for  modal  analysis  and  plotting.  Some  of  its  other  functions 
include  calls  to  various  subroutines  as  follows: 

EIPROP  -  Generates  geometric  data  for  VLM  panels,  and  nodal 
mass  and  flexural  rigidity  distribution  for  finite 
element  analysis. 

VLM  -  Computes  influence  coefficients  for  flat,  swept  wing 
modeled  with  Vortex  Lattice  panels. 

FUNC  -  Sets  right  hand  sides  of  the  dynamical  equations.  If 
time,  T,  equals  zero,  they  are  set  for  initial, 
steady  state  aerodynamic  loads. 

BEAMK  -  Computes  the  reduced  stiffness  matrix  for  a  beam. 

INIT  -  Solves  for  the  initial  wing  deformation 

STRESS  -  Computes  moment  and  shear  at  wing  root. 

DATOUT  -  Displays  results  of  initial  computation  such  as 
flying  speed,  initial  wing  deformation,  and 


NADC  88014-60 


Computer  Program  Descript  ion 


aerodynamic  force  dietribution. 

RK2  -  Advances  the  computation  through  one  full  time-step 
using  the  Runge-Kutta  scheme. 

llie  major  variables  in  NANDYN  include  the  following; 

WO  -  Wing  lateral  deformation  with  respect  to  C.G. 

*1  -  Wing  lateral  velocity  with  respect  to  C.G. 

W2  -  Wing  lateral  acceleration  with  respect  to  C.G. 

BM  -  Nodal  mass 

BC  -  Nodal  damping  (set  to  zero) 

F  -  Right  hand  side  of  dynamical  equation 
BR  -  Reduced  stiffness  matrix 

The  computational  sequence  in  MANDYN  is  as  follows: 

For  unit  air  speed,  the  aerodynamic  force  distribution 
is  established.  For  the  gross  weight  of  the  aircraft, 
required  airspeed  is  determined  and  the  aerodynamic 
forces  are  computed.  l%e  stiffness  matrix  is 
assembled. 

The  steady  state  wing  deformations  are  then  obtained 
and  applied  as  initial  condition  for  WO.  The  vectors 
W1  and  W2  are  set  to  zero.  The  wing  root  shear  and 
moment  are  saved  for  computing  the  dynamic-to-static 
stress  ratio. 

Files  are  then  initialzed  and  data  needed  for  modal 
analysis  are  written.  The  time  step  loop  is  entered 
and  the  displacement,  velocity  and  acceleration 
vectors  are  calculated  by  the  Runge-Kutta  algorithm. 
These  are  recorded  for  post-processing  along  with 
the  dynamic-to-static  stress  ratio. 

The  saved  results  from  MANDYN  in  the  file  RESULT. FIL  can  be 
plotted  after  processing  by  the  Fortran  program  PROCES.  Actual 
plotting  is  done  by  the  Basic  program  MDPLOT;  IBM  graphics  is 
required. 
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The  results  saved  in  EIG.FIL  can  be  processed  by  the  Fortran 
program  EIGVAL  to  produce  natural  frequencies,  mode  shapes  and 
participation  factors.  Weighting  function  for  the  participation 
factors  is  the  steady  state  aerodynamic  force  distribution.  By 
calling  the  Fortran  subroutine  GRAPH,  the  results  are  processed 
for  plotting  by  the  Basic  program  BODES;  IBM  graphics  is  again 
required. 

A  few  words  about  computation  of  the  aerodynamic  loads:  In  the 
subroutine  VLM,  the  influence  coefficients  are  stored  in  the 
augmented  part  of  the  square  matrix  in  the  two  dimensional 
array  W.  These  are  obtained  by  the  subroutine  GAUSS  and 
constitute  the  inverses  of  the  original  VLB  matrix.  At  each 
time  step,  in  the  subroutine  AERODYN,  the  current  wing 
vibrational  velocity  is  used  to  compute  the  vortex  strengths 
and  the  resulting  aerodynamic  forces.  Also,  in  AERODYN  the 
force  vector  F  is  set  to  reflect  the  fuselage  area,  in  an 
arbitrary  way,  to  achieve  the  desired  tuning  for  the  model.  If 
other  cases  are  to  be  analyzed,  this  will  need  to  be  modified. 

In  the  routine  STRESS,  shear  and  moment  can  be  computed  either 
from  the  element  matrix  or  by  summation  of  forces;  the  former 
is  commented  out  of  the  program  because  the  latter  is  more 
efficient;  results  are  the  same  with  either.  Also,  combining 
stresses  by  taking  the  square  root  of  the  sum  of  squares  of 
shear  and  bending  components  is  included  in  MANDYN. 

The  next  two  pages  show  the  user  interaction  with  the  program. 
Input  data  may  be  related  to  various  files  generated  during  a 
run  through  the  Run  No.  which  is  recorded  on  each  file. 
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maneuver  dynamics  RR06RAM  *•« 


Gross  Ut  (lbs)  ~  9Sfi48.000 

wing  group  S  wing  futi  ut  <IDs)  - 
Asptet  Ratio  =  7  270 

Span  (ft):  84  IDO 

wing  root  tnieknttt,  (in  )  :  17  100 

Anglt  of  attack  («ag>  s  5  000 
First  oanding  fraguancy  of  wing  (Hi)  = 
Ratsp-up  tiaa  (Soc)  =  000000 

Load  Factor  (g)  s  1  SO 


5799  000 


4  2S0 


To  cnanga  tna  abova,  kay-in  tna  propar  nu«bar 
fro*  1  thru  8  RETURN  for  no  cnanga 


Piaasa  kay-m  tna  vatua. 


Gross  ut.  (lbs)  :  S8S48.000 

uing  group  t  wing  fual  ut.  (lbs)  : 

Aspact  Ratio  s  7  270 

Span  <f t  )  s  84. 100 

uing  root  tnieknass,  (in.)  *  17  800 

Angla  of  attack  (dag)  s  S  000 

First  banding  fraguaney  of  wing  (Hi)  s 

Ra*p-up  tlMt  (SdC)  :  000000 

Load  Factor  (g)  :  2.  00 


5799  000 


4  280 


To  cnanga  tna  abova,  kay-in  tna  propar  nu»oar 
fro*  I  tnru  8  RETURN  for  no  cnanga 
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9.  RESULTS 


Tlie  aircraft  selected  for  aodeling  is  a  real  life,  variable- 
geometry  airplane.  Hie  configuration  chosen  for  analysis  has 
the  following  parameters: 


Geometry 


Span 

Aspect  Ratio 
Tip  Chord 
Pivot  Chord 

Root  Chord  (Theoretical) 
Tip  Thickness 
Root  Thickness 
Sweep  (LAE) 


64.1  ft. 
7.27 

44.2  in. 
131.5  in. 
167.2  in. 

3.98  in. 
17.8  in. 
20.0  deg. 


Weights 


Gross  wt. 
Wing  group 


58648  lbs. 
5799  lbs. 


1st  Wing  Bending 
2nd  Wing  Bending 
3rd  Wing  Bending 


4.26  hz. 
13.89  hz. 
27.23  hz. 


~  6  in. 


Invoking  symmetry,  only  half  of  the  wing  was  modeled  with  eight 
beam  elements  and  eight  VLM  panels.  The  resulting  nodal 
spacing  was  such  that  a  node  coincided  with  the  pivot  of  the 
variable  geometry  wing.  The  lift  developed  by  each  panel  was 
apportioned  equally  to  the  adjoining  structural  nodes. 
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The  model  was  then  tuned  iteratively  to  reproduce  level  flight 
tip  deflection  of  6  in.  as  well  as  the  three  synnetric  bending 
frequencies.  The  parameters  for  accomplishing  this  were  the 
mass  and  flexural  rigidity  distributions,  and  to  a  lesser 
extent,  the  lift  distribution,  particularly  between  the  pivot 
and  the  wing  root  where  the  fuselage  planform  adds  to  the  wing 
area. 

It  is  appropriate  to  comment  here  on  the  general  problem  of 
synthesizing  simpler  dynamical  models  to  mimic  the  behaviour 
of  more  complex  systems.  In  applied  mechanics  literature, 
this  is  known  as  the  "inverse  problem."  Gladwell  [28]  recently 
reviewed  the  literature  relating  to  the  "unique  reconstruction 
of  a  vibrating  system  from  natural  frequency  data.  The 
literature  is  divided  into  two  groups  -  those  papers 
concerning  discrete  systems,  for  which  the  inverse  problems 
are  closely  related  to  matrix  inverse  eigenvalue  problems,  and 
those  concerning  continuous  systems  governed  either  by  one  or 
the  other  of  the  Sturm-Liouville  equations  or  by  the  Euler- 
Bernoulli  equation  for  flexural  vibrations  of  a  thin  beam." 

An  examination  of  the  reviewed  methodologies,  however, 
disclosed  that  they  could  not  be  readily  applied  to  the 
present  problem  involving  a  "free-free"  boundary  condition. 
Gladwell ’s  concluding  remarks  about  the  state-of-the-art  are 
noteworthy:  "Definitive  answers  to  the  question  of  existence, 
uniqueness,  and  reconstruction  have  been  given  only  for  a 
comparatively  small  class  of  vibrating  systems,  namely  certain 
in-line  discrete  systems  and  one-dimensional  continuous 
systems.  No  significant  results  have  been  obtained  to  date  for 
two-dimensional  systems.  In  the  author’s  opinion  these  results 
will  come  only  when  more  is  known  about  the  qualitative  nature 
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of  ouch  systems,  eg.  how  various  frequencies  change  with 
changes  in  boundary  conditions,  how  eigenmodes  corresponding  to 
different  boundary  conditions  may  be  categorized,  and  what 
eigenfrequency  and  mode  data  are  necessary  and  sufficient  for 
reconstruction.”  Thus,  the  model  synthesis  and  tuning  for  the 
present  problem  was  accomplished  by  intuitive  judgements  based 
on  results  from  iterative  modifications  rather  than  any  formal 
procedure. 

To  tune  the  model,  the  lift  distribution  for  unit  velocity  was 
computed  and  the  aircraft  airspeed  was  then  adjusted  to 
equilibrate  the  total  aerodynamic  force  with  the  gross  weight. 
The  wing  deformations  for  steady  level  flight  were  then 
computed  from  a  static  finite  element  analysis  of  the  wing  for 
the  computed  aerodynamic  load  distribution  with  a  maximum 
flexural  rigidity  of  one.  The  required  maximum  flexural 
rigidity  was  then  adjusted  to  produce  a  tip  deflection  of  about 
6  in.  An  eigenvalue  analysis  of  the  resulting  dynamical  system 
was  performed  to  obtain  the  mode  shapes  and  frequencies.  If  a 
frequency  was  not  satisfactory,  then  by  examining  the  mode 
shapes  intuitive  judgements  were  made  about  how  the  mass  and 
flexural  rigidity  distributions  had  to  be  modified  to  shift  the 
frequencies  towards  the  desired  values.  By  repeating  this 
procedure  several  times,  the  final  model  was  fine  tuned  to  a 
point  where  the  desired  tip  deflection  and  the  first  three 
bending  frequencies  were  matched  excellently. 

An  interesting  finding  of  the  tuning  procedure  was  that  in 
order  to  match  the  third  bending  frequency,  the  location  of 
the  maximum  stiffness  in  the  wing  had  to  be  shifted  from  the 
initially  guessed  central  location  to  the  pivot  area.  Pre¬ 
sumably,  the  pivot  area  is  a  hard  point  in  the  wing  structure. 

-41- 
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Tlie  final  mass  and  flexural  rigidity  distributions  for  the 
model  are  shown  in  figures  9.1  and  9.2,  respectively.  The 
aerodynamic  load  distribution  is  shown  in  fig.  9.3.  Wing 
deformation  for  IG  steady  level  flight  is  shown  in  fig.  9.4. 
(At  the  central  node,  full  values  of  mass  and  lift  forces  are 
shown  in  the  figures.  In  the  analysis,  however,  because  of 
symmetry  boundary  conditions  at  this  node,  only  half  the  value 
of  each  is  used.) 

The  mode  shapes  for  the  first  three  bending  frequencies  are 
illustrated  in  fig.  9.5.  The  modal  participation  factors 
(based  on  the  aerodynamic  force  distribution  for  steady  level 
flight)  are  1186,  112  and  393  for  the  first  to  the  third  modes. 
Thus,  the  third  bending  mode  does  contribute  to  the  response 
moderately. 

Results  from  a  typical  run,  with  a  step-change  in  loading  from 
IG  steady  level  flight  to  2G,  are  illustrated  in  the  subsequent 
figures.  Fig  9.6  shows  the  wing  tip  deflection  as  a  function 
of  time.  Approximately  four  cycles  of  the  1st  bending  mode  at 
a  frequency  of  4.27  hz.  can  be  seen. 

The  wing  tip  velocity  is  shown  in  fig.  9.7.  In  addition  to 
the  first  bending  mode,  contribution  from  the  third,  at  a 
frequency  of  27.75  hz.,  is  apparent  in  the  figure.  The  wing 
tip  acceleration  is  shown  in  fig.  9.8.  Contribution  from  the 
third  mode  is  seen  to  be  quite  pronounced. 

The  displacement,  velocity  and  acceleration  histories  of  the 
wing  root  are  shown  in  figs.  9.9  thru  9.11.  As  is  to  be 
expected,  there  is  a  180  degree  phase  difference  with  the  wing 
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tip  histories.  The  spectral  content  is  similar  to  that  of  the 
wing  tip.  Corresponding  histories  for  the  pivot  are  shown  in 
figs.  9.12  thru  9.14. 

The  non-dimensional  wing  root  dynamic  stress,  obtained  as  the 
ratio  of  dynamic  bending  moment  to  that  of  the  static  moment 
for  IG  loading,  is  plotted  in  fig.  9.15.  The  third  mode  is 
seen  to  contribute  to  the  peak  stress  moderately. 

The  C.G.  G-level  history  is  shown  in  fig.  9.16.  Comparing  the 
peaks  with  those  for  the  stress  ratio  in  fig.  9.15,  a  phase 
difference  greater  than  180  degrees  may  be  observed.  This  is 
in  general  agreement  with  test  data  [30].  Also,  while  the 
G-level  varies  between  a  minimum  of  1.942g  and  a  maximum  of 
2.046g,  the  stress  ratio  peaks  at  2.809. 

Table  9.1  summarizes  the  results  from  ninety  different  runs  in 
which  the  load  factors  and  ramping  times  were  varied.  The 
ratios  of  the  stress  overshoot  to  that  corresponding  to  the 
load  factor  are  tabulated  in  percent. 

Fig  9.17  presents  the  results  in  Table  9.1  in  a  graphical  form. 
The  plots  are  very  reminiscent  of  the  results  for  a  damped 
single  degree  of  freedom  system  subjected  to  similar  ramp 
loading.  (The  dips  in  the  curves  should  occur  at  values  of 
1,2,3,...  for  t/T.  The  particular  plotting  routine  uses  only 
straight  line  segments  to  connect  the  points;  hence,  the 
aberration. ) 

From  recorded  flight  data  for  this  aircraft,  we  know  that  the 
load  factor  N2  can  routinely  change  at  a  rate  of  about  3g/sec. 
during  maneuvers.  A  specific  example  is  as  follows;  Within  0.3 
seconds,  changed  from  1.89  to  2.83  [30].  From  Table  9.1, 
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for  a  load  application  time  of  .3  sec.  and  transition  from  Ig 
to  2g,  the  overstress  is  6.26%.  Recognizing  that  our  model  is 
at  sea  level  and  the  actual  flight  is  at  10,000  ft.,  we  can 
surmise  that  the  actual  overstress  will  be  somewhat  higher 
because  of  the  reduced  aerodynamic  damping  of  wing  vibration 
with  increasing  altitude. 


Fig.  9.3  Aerodynaoic  Load  Distribution 
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Table 


Fig.  9.17  Percent  Stress  Overshoot  Vs. 


Non-dimensional  Ramp-up  Time 
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10.  0(WCLUSI(»iS  AND  REO(»MENDATIC»IS 


Conclusions 

The  model  synthesized  during  the  course  of  this  work  suggests 
strongly  that  dynamic  overstressing,  which  is  unaccounted  for 
in  the  present  design  and  analytical  practices  based  on  the  use 
of  load  factors,  can  be  significant  under  certain  conditions. 

The  significance  arises  primarily  from  the  observation  that  a 
few  percent  increase  in  peak  stress  can  produce  a  much  larger 
reduction  in  fatigue  life. 

Whether  the  load  factor  approach  is  valid  or  not  is  dependent 
on  the  lowest  natural  frequency  of  the  vibratory  modes  brought 
into  play  by  the  particular  set  of  loads.  Also,  the  non- 
dimensional  ramp-up  time  (ratio  of  the  actual  ramp-up  time  to 
the  period  of  the  lowest  mode)  is  an  important  parameter.  If  it 
is  less  than  three,  the  load  factor  approach  can  be  problema¬ 
tical. 

In  assessing  the  effect  of  one  maneuver  on  the  next,  while 
specific  cases  have  not  been  studied,  it  is  nevertheless 
reasonable  to  conclude  that  the  non-dimensional  ramp-up  time, 
based  on  the  time  for  transition  from  one  maneuver  to  the  next, 
will  apply. 

Release  of  stores  is  potentially  the  most  serious  operation 
that  can  be  considered  here.  Because  of  the  much  higher 
overstresses  which  result  from  step  changes  in  loading,  even 
small  stores,  when  released,  can  have  a  significant  impact  on 
fatigue  life. 
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Conclusions  and  Recommendations 


The  C.G.  G-level  is  found  to  be  a  poor  indicator  of  the  non- 
dimensional  stress  ratio,  for  both  phasing  and  magnitude.  This 
raises  a  question  on  how  effective  instrumentation  is,  when 
located  at  the  C.G. ,  for  fatigue  life  monitoring  of  various 
critical  areas. 


The  peak  overstress  vs.  ramp-up  time  data  given  in  this  report 
is  applicable  to  sea  level  conditions.  At  higher  altitudes, 
since  aerodynamic  damping  on  wing  oscillations  is  less,  the 
peak  overstress  will  be  greater.  Hence,  computation  of  this 
data  for  several  altitudes  is  recommended. 

Internal  damping  in  the  structure  is  a  mitigating  factor. 
Therefore,  computation  of  the  data  for  several  damping  values 
is  also  recommended. 

Inclusion  of  fully  unsteady  aerodynamics  will  be  useful. 
However,  simulating  such  effects  of  a  full  scale  aircraft  in  a 
simple  model  is  likely  to  be  difficult.  A  study  of  simple 
models  to  accomplish  this  may  prove  to  be  of  long  term 
benefit. 

A  key  finding  of  the  present  study  is  that  the  non-dimensional 
ramp-up  time  and  the  period  of  the  natural  mode  appropriate  to 
the  loading  are  important  parameters  in  limiting  the  dynamic 
stress  overshoot.  How  best  can  the  present  specifications  be 
supplemented  or  modified  is  an  important  question  that 
warrants  an  in-depth  investigation.  Such  a  study  is 
recommended  to  be  of  considerable  benefit  to  future  military 
aviation. 
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$DEBUG 


!!!!!!!!!!!!!!!!!!!!!!'''*>>*'  I  I  '■■  I  I  I  I  I  I  I  I  I  t  I  I  I  I  I  I  I  I  I  I  I  I 


***  FIXED  WING  AIRCRAFT  MANEUVER  DYNAMIC  ANALYSIS  *** 


Version  2.0 


2-17-87 


Prepared  for  the  NAVAL  AIR  DEVELOPMENT  CENTER 


CONRAD  TECHNOLOGIES,  INC. 


M.Mahadeva  Reddi 


PROGRAM  MANDYN 

COMMON  DT,HAFDT,HAFDT2,RCPDT,T,W0(25) ,W1(25) ,W2(25) ,V(25) 
COMMON/STRUCT/BM(25) ,BC(25) , AK(4 , 4) , BK(25, 25) ,CK(50,50) 
COMMON/AERO/NPAN,Cl,C2,AR,ALF,B,CRD,U,G,F(25) ,W(24,48) 
COMMON/PROP/FLEX(48) , ELEN (48) ,CHORD(48) ,THICK(48) , XI (48 ) , Y1 (48) , 
X2(48) ,Y2(48) ,XM(48) , YM(48) ,ARM(48) 

COMMON/F14/FUSMAS , WINGMS , ANG , CR , CT , ROOTTH , TIPTH 
COMMON/RIGI D/AM , G1 , FR ( 2  5 ) 

COMMON/RAMP/RMPTIM, GLEVEL 
DIMENSION  RK(25) 

DATA  RHO  /. 002378/,  PI  /3. 1415926/ 

DATA  GRWT  /58648./ 

DATA  WINGWT  /5799./ 

DATA  ARATIO  /7.27/ 

DATA  SPAN  /64.1/ 

DATA  ROOT  /17.8/ 

DATA  ALPHA  /5./ 

DATA  OMEGA  /4.26/ 

U=l. 

G*32.2 

ANG*20. 

CR=13.9 
CT=3 . 6 
ROOTTH=17.8 
TIPTH=3.98 
RMPTIM=0 . 

GLEVEL=1.5 
WRITE (*,1) 

FORMAT (IHl,/// 

lOX,'***  MANEUVER  DYNAMICS  PROGRAM  ***'//) 

FORMAT (110) 

FORMAT (GIO. 5) 

WRITE (*,25) GRWT , WINGWT , ARATIO , SPAN , ROOT , ALPHA , OMEGA , RMPTIM , 
GLEVEL 

FORMAT (/ IX, '1.  Gross  Wt.  (lbs)  *',F12.3/ 

IX, '2.  Wing  group  &  wing  fuel  Wt.  (lbs)  =',F12.3/ 

IX, '3.  Aspect  Ratio  ,F9.3/1X, '4.  Span  (ft.)  =',F9.3/ 

IX, '5.  Wing  root  thickness,  (in.)  -',F9.3/ 

IX, '6.  Angle  of  attack  (deg)  *',F9.3/ 

IX, '7.  First  bending  frequency  of  wing  (Hz)  =',F9.2/ 

IX, '8.  Ramp-up  time  (Sec)  =',F12.6/ 

IX, '9.  Load  Factor  (g)  =* ,F10,2/) 

WRITE(*,30) 

FORMAT (/IX, 'To  change  the  above,  key-in  the  proper  number'/lx. 
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'from  1  thru  8.  RETURN  for  no  change.  V) 

READ(*,13)M 

IF(M  .EQ.  0)GOTO  40 

WRITE (*,32) 

F0RMAT(/1X,  'Please  key-in  the  value.'/) 

READ(*,15)AI 
1F(M  .EQ.  1)GRWT=AI 
IF(M  .EQ.  2)WINGWT=AI 
IF(M  .EQ.  3)ARATIO=AI 
IF(M.  EQ.  4)SPAN=AI 
IF(M.  EQ.  5)ROOT=AI 
IF(M.  EQ.  6)ALPHA=AI 
IF(M.  EQ.  7)OMEGA=AI 
IF(M.  EQ.  8)RMPTIM=AI 
IF(M.  EQ.  9)GLEVEIr=AI 
GOTO  20 
AR=ARATIO 
ALF=ALPHA*PI/180. 

B=SPAN 
WRITE (*,42) 

FORMAT (/IX, 'Please  key-in  the  tuning  coefficient  for  El.'/) 

READ(*,15)COEFF 

NEL=8 

NPAN=NEL 

NODES=NEL+l 

NDl=NODES+l 

FUSMAS=0 . 5* (GRWT-WINGWT) /G 
WINGMS=0 . 5*WINGWT/G 
PIMEGA=2 . *PI*OMEGA 
HAFSPN=0.5*SPAN 

EI=FUSMAS* (HAFSPN**3 ) * (PIMEGA**2 ) / (6 . * ( 1 . +GRWT/WINGWT) ) 

EI=EI*COEFF 

AM=WINGMS+FUSMAS 

CALL  ELPROP(NEL) 

Cl=4.*PI*RHO*U 
C2=-U*ALF 
CALL  VLM(NEL) 

T=0. 

DO  44  1=1, NDl 
BC(I)=0. 

V(I)=0. 

W0(I)=0. 

W1(I)=0. 

W2(I)«0. 

CALL  FUNC(ND1,T,W0,W1,W2) 

U-SQRT ( 0 . 5*GRWT/F (NDl ) ) 

Cl=4.*PI*RHO*U 

C2=-U*ALF 

UMPH=U*3600./5280. 

WRITE (*, 61 )UMPH 

FORMAT (/IX, 'Computed  air  speed  (MPH)  =',F10.1/) 
DT=1./(100.*OMEGA) 

WRITE (*, 62) DT 

FORMAT (/IX, 'Computed  time  increment,  DT  =',F10.6/ 

'  To  change  DT,  key-in  1.  Otherwise  key-in  0'/) 


(MPH)  =',F10.1/) 


time  increment,  DT  =',F10.6/ 
key-in  1.  Otherwise  key-in  0'/) 
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READ(*,13)NN 
IF(NN.EQ.O)GOTO  67 
WRITE(*,64) 

64  FORMAT (/IX, 'Please  key-in  the  time- increment'/) 

READ(*,15)DT 

WRITE (*, 65) DT 

65  FORMAT (/IX, 'New  DT  =',F10.6/) 

67  HAFDT=0.5*DT 

HAFDT2=HAFDT*DT 

rcpdt=i./dt 

CALL  BEAMK(NEL,EI) 

WRITE(*,69) 

69  F0RMAT(/1X, 'Please  key-in  the  total  no.  of  time  steps.'/) 

READ ( * , 13 ) NTOT 

WRITE (*,70) 

70  FORMAT (/IX, 'Please  key-in  1  for  terminal  display  of  results,'/ 
1  lx, 'OR,  2  for  saving  on  "RESULT. FI L”.  3  will  give  both.'/) 

READ(*,13)IDPLAY 
IF(IDPLAY  .EQ.  l)GOTO  80 
WRITE(*,72) 

72  FORMAT (/IX, 'Please  key-in  the  step  increment  for  saving.'/) 
READ(*,13)NSAVE 

WRITE (*,73) 

73  FORMAT (/IX, 'Please  key-in  a  Run  No.  to  identify  results'/ 

1  IX, 'stored  in  "RESULT. FIL"  on  the  default  drive.'/) 

READ ( * , 13 ) NRUN 
80  CONTINUE 

85  FORMAT (/IX, 'Thanks!  We  are  flying  now!'///) 

DO  100  1=1, NODES 
100  FR(I)=SQRT(BK(I,I)*BM(I) ) 

CALL  INIT (NODES) 

CALL  FUNC(ND1,T,W0,W1,W2) 

CALL  STRESS ( NODES, S I, SHRI) 

STATIC=SQRT(SI**2+SHRI**2) 

CALL  DATOUT(EI) 

WRITE(*,85) 

IF(IDPLAY  .EQ.  1)G0T0  120 

OPEN ( 4 , FILE= ' RESULT .FIL', FORM= ' UNFORMATTED ' , STATUS= ' OLD ' ) 

OPEN ( 7 , FI LE= ' EIG . FI L ' , FORM= ' UNFORMATTED ' , STATUS= ' OLD ' ) 

WRITE (7) NRUN, NODES, (F(I) ,BM(I) , (BK(I,J) ,J=1, NODES) ,1=1, NODES) 
WRITE ( 4 ) NRUN , NEL , NDl , NTOT 
120  DO  1000  J=l,NTOT 
CALL  RK2(ND1) 

IF(IDPLAY.NE.l)GOTO  150 
WRITE(*,140) J,T 
WRITE(*,126) (W0(I) ,1=1, NODES) 

WRITE(*,126) (W1(I) ,1=1, NODES) 

WRITE(*, 126) (W2 (I) ,1=1, NODES) 

126  FORMAT (5E15. 5) 

140  FORMAT(I7,F15.6) 

GOTO  170 

150  NN=MOD(J,NSAVE) 

IF(NN  .NE.  0)GOTO  170 
CALL  STRESS (NODES, S,SHR) 

DYNAMC=SQRT(S**2+SHR**2) 


3 


NADC  88014^0 


SR=DYNAMC/STATIC 
G2=(32.2+W2(ND1) )/32.2 

WRITE(4) J,T,G2,SR, (WO (I) ,W1 (I) ,W2 (I) ,I=1,ND1) 

IF(IDPLAY  .NE.  3) GOTO  170 
WRITE (*, 140) J,T 
WRITE(*,126) (W0(I) ,1=1, NODES) 

WRITE(*,126) (W1(I) ,1=1, NODES) 

WRITE(*,126) (W2(I) ,1=1, NODES) 

170  CONTINUE 
1000  CONTINUE 

CLOSE ( 4 , STATUS= ' KEEP ' ) 

CLOSE ( 7 , STATUS= ' KEEP ' ) 

STOP 

END 

$DEBUG 

SUBROUTINE  VLM(NEL) 

C 

C  VORTEX-LATTICE  SOLUTION  FOR  FLAT,  SWEPT  WING 

C 

COMMON/ PROP/ FLEX (48) , ELEN (48) , CHORD (48) , THICK (48) , XI (48) , Y1 (48) , 
*  X2 (48) ,Y2 (48) ,XM(48) ,YM(48) 

COMMON/ AERO/NPAN, Cl, C2,AR,ALF,B, CRD, U,G,F( 25) ,W(24,48) 

C 

C  AR  -  ASPECT  RATIO 

C  B  -  SPAN 

C  C  -  CHORD 

C  XM,YM  -  COORDINATES  OF  CONTROL  POINT,  3/4  CHORD 

C  X1,Y1  -  LEFT  COORDINATES  OF  BOUND  VORTEX,  1/4  CHORD 

C  X2,Y2  -  RIGHT  COORDINATES  OF  BOUND  VORTEX,  1/4  CHORD 

C 
C 

C  ASSEMBLE  DOWNWASH  MATRIX 

C 

DO  30  M=1,NEL 
DO  20  N=1,NEL 

W(M,N)=WMN(XM(M) ,YM(M) ,X1(N) ,Y1(N) ,X2(N) ,Y2(N) ) 

K=2*NEL-N+1 

W(M,N)=W(M,N)+WMN(XM(M) ,YM(M) ,X1(K) ,Y1(K) ,X2 (K) ,Y2 (K) ) 

20  CONTINUE 
30  CONTINUE 
C 

C  SET  RIGHT  HAND  SIDES 

C 

DO  50  1=1, NEL 
DO  40  J=1,NEL 
40  W(I, J+NEL)=0. 

50  W(I,I+NEL)=1. 

CALL  GAUSS (W, NEL, NEL) 

RETURN 

END 

C 

C !!!!!!!!!!!!!!!!!!!!!!!  i  !!!!!!!!!!!!!  i  2  !!  i  !!!!  2  !!!!!!!!!!!  2  !  1  !!  2  1  !  1  !!  i  ! 
C 

FUNCTION  WMN(XM,YM,X1,Y1,X2,Y2) 

C 
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110 


COMPUTES  (4PI/GAMMA(N) )*W(M,N) 

WHERE  W(M,N)  IS  VELOCITY  AT  THE  CONTROL  POINT  OF 
THE  M-TH  PANEL  DUE  TO  THE  3  VORTICES  OF  THE  N-TH 


XMNl-XM-Xl 

XMN2=XM-X2 

YMN1=YM-Y1 

YMN2=YM-Y2 

XBV*X2-X1 

YBV=Y2-Y1 

D1=SQRT(XMN1**2+YMN1**2) 

D2 =SQRT ( XMN2 * * 2 + YMN2 * * 2 ) 

CON=l ./ (XMN1*YMN2-XMN2*YMN1) 
C1=(XBV*XMN1+YBV*YMN1)/D1 
C2» (XBV*XMN2+YBV*YMN2 ) /D2 

WMN=CON* ( C1-C2 ) - ( 1 . +XMN1/D1) / YMN1+ ( 1 . +XMN2/D2 ) /YMN2 

RETURN 

END 


SUBROUTINE  GAUSS (A,NEQNS,NRHS) 


SOLUTION  OF  LINEAR  ALGEBRAIC  SYSTEM  BY 
GAUSS  ELIMINATION  WITH  PARTIAL  PIVOTING 
FROM  "ANON."  LIBRARY 

[A]  =  COEFFICIENT  MATRIX 

NEQNS  =  NUMBER  OF  EQUATIONS 

NRHS  =  NUMBER  OF  RIGHT  HAND  SIDES 


RIGHT-HAND  SIDES  AND  SOLUTIONS  STORED  IN 
COLUMNS  NEQNS+1  THRU  NEQNS+NRHS  OF  [A] 

DIMENSION  A(24,48) 

NP=NEQNS+1 

NTOT  =NEQNS+NRHS 

GAUSS  REDUCTION 

DO  150  1=2, NEQNS 

SEARCH  FOR  LARGEST  ENTRY  IN  (I-l)TH  COLUMN 
ON  OR  BELOW  MAIN  DIAGONAL 


IM-I-1 

IMAX=IM 

AMAX=ABS(A(IM,IN) ) 

DO  110  J=I, NEQNS 

IF(AMAX.GE.ABS(A(J,IM)))GOTO  110 
IMAX=J 

AMAX=ABS(A(J,IM) ) 

CONTINUE 


SWITCH  (I-l)TH  AND  IMAXTH  EQUATIONS 
IF(IMAX.NE.IM)GOTO  140 
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130 

C 

c 

c 

c 

140 


150 

C 

C 

C 


200 

210 

220 


C 

C! J  1  I  J  I 

C 

C 


* 


c 


20 


1 

2 

3 

4 

5 


30 


1 


50 

55 


DO  130  J=IM,NT0T 

TEMP*A(IM,J) 

A(IM,J)»A(IMAX,J) 

A(IMAX,J)«TEMP 

CONTINUE 


ELIMINATE  (I-l)TH  UNKNOWN  FROM 
ITH  THRU  (NEQNS)TH  EQUATIONS 

DO  150  J=I,NEQNS 
R=A(J,IM)/A(IM,IM) 

DO  150  K=I,NTOT 
A(J,K)=A(J,K)-R*A(IM,K) 

BACK  SUBSTITUTION 

DO  220  K»NP,NTOT 

A ( NEQNS , K ) «A ( NEQNS , K) / A ( NEQNS , NEQNS ) 

DO  210  L=2, NEQNS 

I=NEQNS+1-L 

IP=I+1 

DO  200  J=IP, NEQNS 
A(I,K)=A(I,K)-A(I,J) *A(J,K) 

A(I,K)=A(I,K)/A(I,I) 

CONTINUE 

RETURN 

END 

I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I »• I  I  M  1  !!!!  1  1  I  !  I  I  J I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I 


SUBROUTINE  DATOUT(EI) 

COMMON  DT,HAFDT,HAFDT2,RCPDT,T,W0(25) , W1 ( 25) , W2 ( 25 ) , V ( 25) 
COMMON/STRUCT/BM(25) ,BC(25) ,AK(4,4) , BK (25 , 25) , CK ( 50 , 50) 
COMMON/AERO/NPAN,Cl,C2,AR,ALF,B,CRD,U,G,F(25) ,W(24,48) 
COMMON/PROP/FLEX(48) ,ELEN(48) ,CHORD(48) ,THICK(48) ,X1(48) ,Y1(48)  , 
X2(48) ,Y2(48) ,XM(48) , YM(48) 

COMMON/F14/FUSMAS , WINGMS , ANG, CR, CT , ROOTTH , TIPTH 
COMMON/RIGI D/ AM , G 1 , FR ( 2 5 ) 

WRITE (*, 2 0 ) ANG , CR , CT , ROOTTH , TIPTH , El 
FORMAT (/'  Sweep  Angle  (deg.)  *',F12.2/ 

'  Root  chord  (ft.)  =',F12.3/ 

’  Tip  chord  (ft.)  ■'',F12.3/ 

*  Root  thickness  (in.)  »',F12.2/ 

'  Tip  thickness  (in.)  »',F12.2/ 

*  Flex,  rigidity  mult.  =',E12.5/) 

WRITE (*,30) (I,X1(I) ,Y1(I) ,X2(I) ,Y2(I) ,XM(I) ,YM(I) ,I=1,NPAN) 
FORMAT (/'  Panel' ,5X, 'XI' ,10X, 'Yl' ,10x, 'X2' ,10X, 'Y2' ,10X, 'Xm' , 
lOX, 'Ym'/(I4,6F12.2) ) 

NODES=NPAN+l 

WRITE(*,50) (F(I) ,1=1, NODES) 

FORMAT(/'  Steady-state  aerodynamic  forces  (lbs. ) : '// (5E15 . 5) ) 
WRITE(*,55) (BM(I) ,1=1, NODES) 

FORMAT (/'  Wing  mass  distribution  (slugs) : '// (5E15 . 5) ) 
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$DEBUG 

C 

C!  !  !  !  !  ! 
C 


WRITE(*,60) (FLEX(I) ,I-1,NPAN) 

FORMAT(/^  Wing  flexural  rigidity  distribution: V/(5E15. 5) ) 
WRITE (*,80) (WO (I) ,1*1, NODES) 

FORMAT(/'  Steady-state  wing  deflections  (ft.) s V/(5E15.5) ) 

RETURN 

END 


!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I !!!!!!!!!!!!!!!!!!!!  M  !!!!!!!!!!!  n 

SUBROUTINE  BEAMK(NEL,EI) 

ASSEMBLES  BEAM  STIFFNESS  MATRIX 

NEL  -  NUMBER  OF  ELEMENTS 
AK  -  EliEMENT  STIFFNESS  MATRIX 

BK  -  REDUCED  BEAM  STIFFNESS  MATRIX 

CK  -  BEAM'S  FULL  STIFFNESS  MATRIX 

COMMON/STRUCT/BM(25) ,BC(25) ,AK(4,4) ,BK(25,25) ,CK(50,50) 
COMMON/PROP/FLEX(48)  ,EU:N(48)  ,CH0RD(48)  ,THICK(48)  ,X1(48)  ,Y1(48)  , 
X2(48) ,Y2(48) ,XM(48) ,YM(48) ,ARM(48) 

ILAST*2*(NEL-l)+4 
DO  10  I*1,ILAST 
DO  10  J*1,I 
CK(I, J)=0. 

ACCUMULATE  LOWER  TRIANGLE 

DO  30  N*1,NEL 

EL=ELEN(N) 

ELEI=EI*FLEX(N) 

CALL  ELSTIF(AK,ELEI,EL) 

IP*2*(N-1) 

DO  20  1=1,4 
K=IP+I 
DO  20  J=1,I 
L*IP+J 

CK(K,L)=CK(K,L)+AK(I,J) 

CONTINUE 

FILL  UPPER  TRIANGLE 

DO  40  I*2,ILAST 
11*1-1 

DO  40  J=1,I1 
CK(J,I)»CK(I,J) 

CONTINUE 

APPLY  SYMMETRY  BOUNDARY  CONDITION 
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DO  46  I«1,I1AST 
CK(I,IIAST)»0. 

46  CK(IIiAST,I)-0. 

CK  ( HAST ,  HAST) -1 . 

ELIMINATE  ROTATORY  DOF^S 

11- HAST-1 

12- HAST-2 

DO  100  1-2,12,2 
DIAG-1./CK(I,I) 

DO  70  11-1,11 
IF(II.EQ.I)GOTO  70 
CON=CK(II,I) 

DO  60  JJ»1,I1 

CK(II,JJ)-CK(II,JJ)-CK(I,JJ)*CON*DIAG 
60  CONTINUE 
70  CONTINUE 
100  CONTINUE 

FILL  REDUCED  MATRIX 

IP-NEL+1 
DO  200  1=1, IP 
11-2*1-1 
DO  200  J-1,IP 
JJ=2*J-1 

BK(I,J)-CK(II,JJ) 

200  CONTINUE 

OM=2.*3.1416*4.23 
PERCNT-.02 
DO  300  1=1, IP 

300  BC(I)=PERCNT*2.*BM(I) *OM 

RETURN 

END 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

SUBROUTINE  ELSTIF (AK, El , EL) 

COMPUTES  STIFFNESS  MATRIX  FOR  BEAM  ELEMENT 

AK  -  ELEMENT  STIFFNESS  MATRIX 
El  -  FLEXURAL  RIGIDITY 
EL  -  ELEMENT  LENGTH 

DIMENSION  AK(4,4) 


FILL  LOWER  TRIANGLE 


C-2.*EI/EL**3 
AK(1,1)-6.*C 
AK(2,1)— 3.*C*EL 
AK(2,2)=2.*C*EL*EL 


V’ 


I 

V, 


z* 


I 


I 

Si’ 

k 

il 

s 
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AK(3,1)— AK(1,1) 

AK(3,2)«-AK(2,1) 

AK(3,3)»AK(1,1) 

AK(4,1>-AK(2,1) 

AK(4,2)»0.5*AK(2,2) 

AK(4,3)»AK(3,2) 

AK(4,4)-AK(2,2) 

RETURN 

END 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
SUBROUTINE  ELPROP(NEL) 

GENERATES  ELEMENT  DATA 

COMMON/PROP/FLEX (48) ,ELEN(48) ,CHORD(48) ,THICK(48) ,X1(48) ,Y1(48)  , 
*  X2(48) ,Y2(48) ,XM(48) ,YM(48) ,ARM(48) 

COMMON/STRUCT/BM(25) ,BC{25) ,AK(4,4) ,BK(25,25) ,CK(50,50) 

COMMON/ AERO/NPAN, Cl, C2,AR,ALF,B, CRD, U,G,F( 25) ,W(24,48) 
COMMON/F14/FUSMAS , WINGMS , ANG , CR , CT , ROOTTH , TIPTH 
DIMENSION  VOL (24) 

DATA  PI  /3. 1415926/ 

RTH-ROOTTH/20. 

TTH-TIPTH/2 . 

BHAF»B/2. 

EL=BHAF/FLOAT (NEL) 

CRMCT»CR-CT 
TANS«-TAN(ANG*PI/180. ) 

Yl(l)— BHAF 
Y2(1)=-BHAF+EL 
YM ( 1 ) *-BHAF+0 , 5*EL 
ARM(1)=BHAF 

NEL2»2*NEL 
DO  10  I=2,NEL2 
Y1(I)-Y1(I-1)+EL 
Y2(I)»Y2(I-1)+EL 
ARM(I)-ARM(I-1)-EL 
10  YM(I)-YM(I-1)+EL 


DO  20  I-1,NEL 

XI (I) -Y1 (I) *TANS+0. 25* (CT+CRMCT* ( Y1 (I) +BHAF) /BHAF) 
X2 (I) «Y2 (I) *TANS+0. 25* (CT+CRMCT* (Y2 (I) +BHAF)/BHAF) 
CHORD ( I ) -CT+CRMCT* ( YM ( I ) +BHAF) /BHAF 
THI CK ( I ) -TTH+ ( RTH-TTH ) * ( YM ( I ) +  BHAF ) / BHAF 
XM ( I ) -YM ( I ) *TANS+0 . 7 5*CHORD ( I ) 

FLEX(I)-1. 

20  ELEN (I) -EL 

NELl-NEL-1 
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DO  30  I«1,NEL1 

DUM-FLOAT ( I - 1 ) / FLOAT ( NELl - 1 ) 

30  FLEX(I)»0.35+0.65*DUM 

C 

NELl-NEL+1 

DO  35  I=NEL1,NEL2 

XI (I) *-Yl (I) *TANS+0. 25* (CT+CRMCT* (Y1 (I) -BHAF)/BHAF) 
X2 (I) »-Y2 (I) *TANS+0. 25* (CT+CRMCT* (Y2 (I) -BHAF)/BHAF) 
CHORD ( I ) -CT+CRMCT* ( YM ( I ) -BHAF) /BHAF 
35  XM(I)— YM(I)*TANS+0.75*CHORD(I) 


C 


v»o. 

DO  40  I-1,NEL 

VOL ( I ) -CHORD ( I ) *ELEN ( I ) *THICK ( I ) 
BM(I)-0. 

40  V=V+VOL(I) 

BM(NEL+l)-0. 


DO  50  1=1, NEL 
EIMASS-0 . 5*VOL ( I ) *WINGMS/V 
BM ( I ) -BM ( I ) +ELMASS 
50  BM(I+l)-BM(I+l)+EIilASS 
FDIV-FUSMAS/15. 

BM (NEL) =BM (NEL) +FDIV*10 . 

BM(NEL+1) -BM(NEL+1) +FDIV*5. 

BM (NEL+2 ) -WINGMS+FUSMAS 

Y1(1)=Y1(1)+0.75*EL 

YM(1)=YM(1)+0.375*EL 

Y2 (NEL2 ) =Y2 (NEL2 ) -0 . 7  5*EL 

YM(NEL2)=YM(MEL2)-0.375*EL 

XI ( 1 ) =Y1 ( 1 ) *TANS+0 . 2  5* ( CT+CRMCT* ( Y1 ( 1 ) +BHAF) /BHAF) 

X2 ( 1 ) -Y2 ( 1 ) *TANS+0 . 2  5  * ( CT+CRMCT* ( Y2 ( 1 ) +BHAF ) /BHAF ) 

CHORD ( 1) -CT+CRMCT* ( YM ( 1) +BHAF) /BHAF 
XM ( 1 ) -YM ( 1 ) *TANS+0 , 7  5*CHORD ( 1 ) 

XI (NEL2 ) — Y1 (NEL2 ) *TANS+0 .25* (CT+CRMCT* ( Y1 (NEL2 ) -BHAF) /BHAF) 
X2 (NEL2 ) — Y2 (NEL2 ) *TANS+0 . 25* (CT+CRMCT* ( Y2 (NEL2 ) -BHAF) /BHAF) 
CHORD (NEL2 ) -CT+CRMCT* ( YM (NEL2 ) -BHAF) /BHAF 
XM (NEL2 ) — YM (NEL2 ) *TANS+0 . 75*CHORD (NEL2 ) 

BM(l)-63./G 

BM(2)=109./G 

BM(3)-177./G 

BM(4)-343./G 

BM(5)-410./G 

BM(6)-457./G 

BM(7)-8038./G 

BM(8)-9648./G 

BM(9)-10079./G 

BM(10)-29324./G 

FLEX(l)-  .1 

FLEX (2)-  .3 

FLEX(3)-  .325 

FLEX(4)-  .35 

FLEX(5)-1. 

FLEX(6)-1. 

FLEX(7)-  .8 
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FLEX(8)«  .8 

FUM* 2 1 0 7 . / ( F LO AT ( NE L - 3 ) * 3 2 . 2 ) 

DO  60  1*4, NEL 
60  BM(I)«BM(I)+FUM 
RETURN 
END 

DEBUG 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!!!!!!!!!!!!!!!!!!!!!!!!!! 

SUBROUTINE  FUNC(N,T,W0,W1,W2) 

DIMENSION  W0(1) ,W1(1) ,W2(1) 

COMMON/STRUCT/BM(25) ,BC(25) ,AK(4,4) ,BK(25,25) ,CK(50,50) 

COMMON/ AERO/NPAN, Cl, C2,AR,ALF,B, CRD, U,G,F( 25) ,W(24,48) 
COMMON/RIGID/AM , G1 , FR ( 2 5 ) 


;? 

'I* 

♦'l 


N1«N-1 


CALL  AERODN(N,T,Wl) 

W2 (N)=F(N)/BM(N)-G 
IF(T.EQ.O. )W2 (N)-0. 

DO  100  1*1, N1 

W2(I)=F(I)-BC(I)*W1(I)-BM(I)*(G+W2(N)  ) 

DO  50  J*1,N1 

50  W2(I)*W2(I)-BK(I,J)*W0(J) 

100  W2(I)=W2(I)/BM(I) 

RETURN 

END 

I  I  I  I  I  I  I  I  I  I  I  I  )  I  I  I  I  I  I  «  I  M  t  t  «  1  I  I  t  t  t  I  I  I  I  I  «  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  t  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I 


SUBROUTINE  AER0DN(N,T,W1) 

COMMON/ PROP/ FLEX (48) , ELEN (48) ,CHORD(48) ,THICK(48) , XI (48 ) , Y1 (48 ) , 
X2(48) ,¥2(48) ,XM(48) ,YM(48) ,ARM(48) 
COMMON/AERO/NPAN,Cl,C2,AR,ALF,B,CRD,U,G,F(25) ,W(24,48) 
COMMON/RIGID/ AM , G1 , FR ( 2  5 ) 

DIMENSION  Wl(l) ,GAMMA(24) ,V(24) ,FA(24) 

COMMON/RAMP/RMPTIM, GLEVEL 

N1*N-1 

FAC*1 . 0 

IF(T.LE. .005)GOTO  45 

FAC*1.+(GLEVEL-1.)*(T-.005)/(RMPTIM-.005) 


45 

50 

100 


IF (RMPTIM . EQ . 0 . ) FAC*GLEVEL 
IF (T . GT . RMPTIM) FAC*GLEVEL 
F(1)*0. 

DO  50  J*2,N1 
F(J)*0. 

V(J-1)-ELEN(J-1)*C1*(FAC*C2+0.5*(W1(J-1)+W1(J) ) ) 

DO  100  I*1,NPAN 

FA(I)*0. 

DO  100  J*1,NPAN 
JN*J+NPAN 

FA(I)*FA(I)+W(I,JN) *V(J) 

F(N)*0. 
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IX}  200  J«1,NPAN 
F(N)-F(N)+FA(J) 

HAF-0.5*FA(J) 

F(J)«F(J)+HAF 
200  F(J+1)-F(J+1)+HAF 

F(10)-F(10)+2.*(F(9)+F(8) )+F(7) 

F(9)«F(9)*3. 

F(8)«F(8)*3. 

F(7)»F(7)*2. 

RETURN 

END 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!i!!!!!!! 
SUBROUTINE  STRESS (NODES, S,SHR) 

COMMON  DT,HAFDT,HAFDT2,RCPDT,T,W0(25) ,W1(25) ,W2(25) ,V(25) 
COMMON/STRUCT/Wl(25) ,BC(25) , AK(4 , 4) , BK(25 , 25) ,CK(50,50) 
COMMON/PROP/FLEX (48) ,ELEN(48) ,CHORD(48) ,THICK(48) ,X1(48) , Yl(48) , 
*  X2(48) ,Y2(48) ,XM(48) ,YM(48) , ARM (48) 

COMMON/ AERO/NPAN, Cl, C2,AR,ALF,B, CRD, U,G,F( 25) ,W(24,48) 

COMMON  /RIGID/AM, G1,FR(25) 

S=0 . 5* (8 . *W0 (NODES-1) -7 . *W0 (NODES) -WO (NODES -2 ) ) 

RETURN 

END 

Shear  and  moment  computation  from  pseudo-forces  follows: 
Nl=NODES-l 
SHR=0 . 

S=0. 

DO  10  1*1, N1 

DUM=F(I)-BM(I) *(W2(I)+G+W2(NODES+l) ) 

SHR-SHR+DUM 
10  S=S+DUM*ARM(I) 

RETURN 

END 

Moment  computation  from  element  follows: 

Nl*NODES-l 
NTHETA»2*N1 
THETA- 0 . 

DO  10  J-1, NODES 
NC»2*J-1 

1 0  THETA-THETA+CK ( NTHETA , NC ) *W0 ( J ) 

THETA— THETA/ CK  (NTHETA ,  NTHETA) 

S-AK (4,1)*(W0(N1)-W0( NODES ) ) +AK (4,2) *THETA 

RETURN 

END 

!!!!!!!!!!!!!!!!!!!!!!!! i i !! 1 !!!!!!!!  1  !  M  !!!!!  M  !!!!!!!!!!!!!!!!!!!!!!  1 
SUBROUTINE  INIT (NODES) 

COMMON  DT,HAFDT,HAFDT2,RCPDT,T,W0(25) ,W1(25) ,W2 (25) ,V(25) 
COMMON/STRUCT/BM(25) ,BC(25) , AK(4 , 4) , BK(25 , 25) ,CK(50,50) 

COMMON/ AERO/NPAN, Cl, C2,AR,ALF,B, CRD, U,G,F( 25) ,W(24,48) 
COMMON/RIGID/ AM , G1 , FR ( 2  5 ) 


. 'afv***  Ia* '*aA'*aA  4.# 


ft  ft  A»ft.«»'t  ■»!  ■*<  »ift'’.ift  ,*ft^Jtft  ,«4*^f«*^iB*-***  ■«*  ftw>tft«ft  ft.*  ^’ftil'i.i  i*«  **-  «WA»»*fttA  ft»A*ft 
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DIMENSION  DK(24,48) 


NDl-NODES+1 

CALL  AER0DN(ND1,T,W1) 


DO  10  I«l, NODES 
DK(I,ND1)-F(I)-BM(I)*G 
DO  10  J«l, NODES 
10  DK(I,J)«BK(I,J) 


CALL  GAUSS ( DK, NODES, 1) 


DO  20  1*1, NODES 

20  WO(I)«DK(I,ND1)-DK(NODES,ND1) 


RETURN 

END 


III  I  I  I t  I  I  I  I  I  I  I  I  I  I  t  I  I  I  I  I  I  t  I  I  I  I  I  I  I  I  I  I  I  I  I 
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SUBROUTINE  RK2 (N) 


RUNGE-KUTTA  FOR  DIFF.  EQS.  OF  2ND  ORDER 


F(T,W,W  ) 


COMMON  DT,HAFDT,HAFDT2,RCPDT,T,W0(25) ,W1(25) ,W2 (25) ,V(25) 
DIMENSION  Y0(25) ,Y1(25) ,Y2(25) 

DIMENSION  AKN(25) ,AK(25) ,AKP(25) 


INITIALLY,  SETUP  THE  FOLLOWING: 


DT,  T,  WO,  W1 

HAFDT*0.5*DT 

HAFDT2-0.5*DT**2 

rcpdt*i./dt 

W2,  BY  —  CALL  FUNC(N,T,W0,W1,W2) 
V*DT*W1 


DO  10  1*1, N 

AKN ( I ) -HAFDT2 *W2 ( I ) 

AK(I)-AKN(I) 

AKP(I)-AKN(I) 


Tl-T+HAFDT 


DO  20  1*1, N 

Y0(I)-W0(I)+0.5*V(I)+O.25*AK(I) 


DO  50  III«1,2 


DO  30  1*1, N 

Y1(I)=(V(I)+AKN(I) )*RCPDT 


■jAOjjVijjVjjjjjjjWjjjWwWgggjjWg 
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CALL  FUNC(N,T1,Y0,Y1,Y2) 


DO  40  I-1,N 

AKN ( I ) -HAFDT2 * Y2 ( I ) 

AK(I)-AK(I)+AKN(I) 

AKP(I)»AKP(I)+2.*AKN(I) 


CONTINUE 


T=T+DT 
DO  60  1*1, N 

Y0(I)-W0(I)+V(I)+AKN(I) 
Y1(I)»(V(I)+2.*AKN(I) )*RCPDT 


CALL  FUNC(N,T,Y0,Y1,Y2) 


DO  70  I-1,N 

W0(I)*W0(I)+V(I)+AK(I)/3. 
V(I)*V(I)+(AKP(I)+HAFDT2*Y2(I) )/3. 
W1(I)*V(I)*RCPDT 


CALL  FUNC(N,T,W0,W1,W2) 


RETURN 

END 


$DEBUG 

I  I  I  I  I  I 
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Processes  results  from  MANDYN  program  for  plotting  by  MDPLOT.BAS 


PROGRAM  PROCES 


DIMENSION  Y(3,27) ,YS(3,27) ,YL(3,27) , 

NY (3, 27) ,YAMP(3,27) ,SCL(3,27) 


WRITE(*,20) 

FORMAT(//'  Please  input  first  and  last  step  numbers  for  plot.'/) 


READ(*,25)NB 
READ(*,25)NL 
FORMAT (110) 


OPEN ( 4 , FILE* ' RESULT . FIL' , FORM* 'UNFORMATTED ' , STATUS* ' OLD ' ) 


READ ( 4 ) NRUN , NEL, NODES , NTOT 


NDl=NODES+l 
DO  30  1=1,3 
DO  30  J=1,ND1 
YS(I,J)=l.e+20 
YL(I,J)*-l.e+20 


DO  60  ISTEP=1,NL 


•ji  ‘k*  'l».  jf.  *».  *i. 
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READ(4) J,y(l,NDl) ,Y(2,ND1) ,y(3,NDl) , 

(y(i,i) ,y(2,i),y(3,i) ,i=i, nodes) 

IF(ISTEP.LT.NB)GOTO  60 

DO  50  11*1,3 
DO  50  JJ»1,ND1 

ys(ii,jj)-AMiNi(ys(ii,jj) ,y(ii,jj) ) 
yL(ii,jj)«AMAxi(yL(ii,jj) ,y(ii,jj)) 

CONTINUE 

CLOSE ( 4 , STATUS- ' KEEP ' ) 

OPEN ( 4 , FILE- ' RESULT . FIL^ , FORM- ^UNFORMATTED ' , STATUS- ' OLD ' ) 
OPEN ( 7 , FILE- ' PLOT . FIL' , FORM- ' FORMATTED ' , STATUS- ' NEW • ) 

READ ( 4 ) NRUN , NEL , NODES , NTOT 

WRITE (7,25) NRUN , NL , NODES 

WRITE (7, 83) ( (yS(I,J) ,J=1,ND1) ,1=1,3) 

WRITE (7, 83) ( (yL(I,J) ,J=1,ND1) ,1=1,3) 

NS AVE- 1+ (NL-NB)/400 
NPOINT=FLOAT (NL-NB) /FLOAT (NSAVE) 

XAMP=639 . /FLOAT (NPOINT) 

DO  70  1=1,3 
DO  70  J=1,ND1 
CC»yL(I,J)-yS(I,J) 

IF(CC.LT.1.E-10)CC=1. 
yAMP(I,J)=180./CC 
SCL(I, J)=CC/180. 

CONTINUE 

WRITE (7, 83) ( (SCL(I,K) ,K=1,ND1) ,1=1,3) 

XX=0. 

DO  90  ISTEP=1,NL 

READ(4) J,y(l,NDl) ,y(2,NDl) ,y(3,NDl) , 

(y(i,i) ,y(2,i),y(3,i) ,i=i, nodes) 

IF(ISTEP.LT.NB)GOTO  90 
NN-MOD ( ISTEP , NSAVE ) 

IF(NN.NE.O)GOTO  90 

XX-XX+XAMP 

NX-XX 

DO  80  1=1,3 
DO  80  K-1,ND1 

Ny(i,K)-200.-(y(i,K)-ys(i,K) )*yAMP(i,K) 

WRITE (7, 25) NX 

WRITE (7, 82) ((Ny(I,K) ,K-1,ND1) ,1=1,3) 


Lh' 
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FORMAT (1615) 

FORMAT (5E13. 4) 

CONTINUE 

CLOSE ( 4 , STATUS= ' KEEP ' ) 
CLOSE ( 7 , STATUS= ' KEEP ' ) 

STOP 

END 


PROGRAM  EIGVAL 

COMPUTES  MODE  SHAPES  AND  FREQUENCIES 

NEEDS  FULL,  SYMMETRIC  STIFFNESS  AND  DIAGONAL  MASS  MATRIC 

DIMENSION  A(30,30) ,B(30) ,U(30,30) ,W(30) ,KLUE(6) , PF{30) ,GM(30) , 
AM(30) ,F(30} 

CHARACTER  TITLE1*12 ,TITLE2*12 
DATA  TITLE 1/'  MODE  SHAPE  V 
DATA  TITLE2/'FREQ.  IN  HZ-V 
DATA  KLUE/2,2,1,1,1,0/ 

OPEN ( 4 , FILE= ' EIG . FILS F0RM= 'UNFORMATTED STATUS= ' OLD ' ) 
READ(4)NRUN,N, (F(I) ,B(I) , (A(I,J) ,J=1,N) ,I=1,N) 

DO  5  1=1, N 
AM(I)=B(I) 

CALL  EIGEN(A,B,U,N,R1) 

DO  15  J=1,N 
PF(J)=0. 

GM(J)=0. 

DO  10  1=1, N 

GM(J)*GM(J)+AM(I)*U(I,J)**2 

PF(J)=PF(J)+U(I,J)*F(I) 

PF(J)*ABS(PF(J) ) 

WRITE (*, 30) R1 
FORMAT (F2 0.0) 

DO  40  J=1,N 
AAA=ABS(B(J) ) 

VALUE=SQRT (AAA) /6. 2831852 
WRITE(*,49)B(J) , VALUE, PF(J) 

WRITE(*,50) (U(I,J) ,I=1,N) 

FORMAT (3F16. 2, ) 

FORMAT (5E15. 6) 

CLOSE ( 4 , STATUS- ' KEEP ' ) 

DO  90  K-2,6 
CON-1. 

IF(U(1,K)  .LT.0.)CON— 1. 

DO  60  1=1, N 
W(I)=CON*U(I,K) 

VALUE=SQRT(B(K) )/6. 283 1852 

CALL  GRAPH (W,N, TITLE 1,TITLE2,NRUN, VALUE, KLUE(K) ) 

STOP 

END 

COMPUTES  EIGENVALUES  BY  THE  CYCLIC  JACOBI  METHOD 


m 
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SUBROUTINE  EIGEN (A, B,U,N,R1) 

DIMENSION  A(30,30) ,B(30) ,U(30,30) 

SI  -  NUMBER  OF  SIGNIFICANT  DIGITS 

Sl=6. 

Z=2.*S1 

T1=1./(10.**Z) 

R=5*N*N 

R1=0. 

T2=0.1 
N1=N-1 
DO  20  1=1, N 
B(I)=1./SQRT(B(I)  ) 

DO  20  J=1,N 
U(I,J)=0. 

DO  30  1=1, N 
U(I,I)=1. 

DO  30  J=1,N 

A(I,J)=B(I)*A(I,J)*B(J) 

PERFORM  ONE  CYCLE  OF  ROTATIONS 

CALL  ROTATE(A,U,Rl,Xl,N,Nl,T2) 

IF  (Xl.LT.Tl)GOTO  50 

IF(Rl.GT.R)GOTO  41 

T2=0.1*X1 

GOTO  40 

WRITE (*, 45) R1 

FORMAT('  ***NO  CONVERGENCE  AFTER' , FI 5.6, '  ROTATIONS***') 

STOP 

CALL  NORM(A,B,U,N,Nl) 

RETURN 

END 

SUBROUTINE  NORM(A, B,U,N,N1) 

DIMENSION  A(30,30) ,B(30) ,U(30,30) 

NORMALIZE  EIGENVECTORS 

DO  70  1=1, N 
DO  70  J=1,N 
U(I,J)-U(I,J)*B(I) 

DO  80  1=1, N 
B(I)»A(I,I) 

ORDER  SOLUTIONS 

DO  200  1=1, N1 
11=1+1 
2=B(I) 

M=I 

DO  110  J=I1,N 
IF(Z.LT.B(J) )GOTO  110 
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Z*B(J) 

M=J 

CONTINUE 

B(M)=B(I) 

B(I)=Z 

DO  150  J=1,N 
Z=U(J,I) 

U(J,I)=U(J,M) 

U(J,M)=Z 

CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUTINE  ROTATE (A,U,R1,X1,N,N1,T2) 
DIMENSION  A(30,30) ,U(30,30) 

X1=0. 

DO  2000  K=1,N1 
K1=K+1 

DO  1000  L=K1,N 
A1=A(K,K) 

A2=A(K,L) 

A3=A(L,L) 

X=A2*A2/(A1*A3) 

CHECK  IF  ROTATION  IS  NECESSARY 

IF(X.GT.Xl)GOTO  90 

GOTO  95 

X1=X 

IF(X.LT.T2)GOTO  1000 
R1=R1+1. 

COMPUTE  ANGLE 

IF(Al.EQ.A3)GOTO  100 

Z=0.5*(A1-A3)/A2 

Z1=1.+1./(Z*Z) 

T=-Z*(1.+SQRT(Z1) ) 

GOTO  110 
T=l. 

C=1./SQRT(1.+T*T) 

S*C*T 

S2=S*S 

C2=C*C 

A(K,L)=0. 

TRANSFORM  DIAGONAL  ELEMENTS 


A0=2.*A2*C*S 

A(K, K) =A1*C2+A0+A3*S2 

A(L,L)=A1*S2-A0+A3*C2 


TRANSFORM  OFF-DIAGONAL  ELEMENTS 
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DO  900  1*1, N 
IF(I.LT.K)GOTO  160 
IF(I.GT.K)GOTO  180 
GOTO  900 
160  A0=A(I,K) 

A ( I , K) =C*A0+S  *A ( I , L) 

A ( I , L) =-S*A0+C*A ( I , L) 

GOTO  900 

180  IF(I.LT.L)GOTO  190 
IF(I.GT.L)GOTO  200 
GOTO  900 
190  A0=A(K,I) 

A(K,I)=C*A0+S*A(I,L) 

A ( I , L) =-S*A0+C*A ( I , L) 

GOTO  900 
200  A0=A(K,I) 

A(K,I)=C*A0+S*A(L,I) 

A(L,I)=-S*A0+C*A(L,I) 

900  CONTINUE 

GENERATE  EIGENVECTORS 

DO  950  1=1, N 
U0=U(I,K) 

U ( I , K) =C*U0+S*U ( I , L) 

950  U(I,L)=-S*U0+C*U(I,L) 

1000  CONTINUE 
2000  CONTINUE 
RETURN 
END 

SUBROUTINE  GRAPH ( W , N , TITLEl , TITLE2 , NRUN , VALUE , NFI L) 

SMOOTH  CURVE  FOR  N  EQUALLY  SPACED  POINTS 
NFIL  =  3  FOR  ONE  PLOT  ONLY 

NFIL  =  2  FOR  FIRST,  WITH  MORE  PLOTS  TO  COME 
NFIL  =  1  FOR  MORE  PLOTS  TO  COME 
NFIL  =  0  FOR  LAST 

DIMENSION  W(30) ,WP(30) ,NY(640) 

CHARACTER  TITLE1*12 , TITLE2*12 
C 

N1=N-1 
DO  10  1=2, N1 

10  WP(I)=0.5*(W(I+1)-W(I-1) ) 

WP(N)=0. 

C  ; COERCE  SYMMETRY  ABOUT  LAST  POIN 

WMAX=-l.E+25 
WMIN=l.E+25 
DO  20  1=1, N 
WMAX=AMAX1 (WMAX , W ( I ) ) 

20  WMIN=AMIN1(WMIN,W(I) ) 

AMP=WMAX-WMIN 
IF(AMP.EQ.O. )AMP=1. 

ORG=100 . *WMIN/AMP 


o  o 
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NORG=150.+ORG 

NX=639/N1 

MX=0 

A1=2.*(W(2)-W(1))-WP(2) 

A2=WP(2)-W(2)+W(1) 

; PARABOLIC  FOR  FIRST  INTERVAL, 
MATCH  SLOPE  AT  2ND  POINT 

DX=1. /float (NX) 

X»-DX 

DO  30  1=1, NX 

X=X+DX 

MX=MX+1 

y=W(l)+Al*X+A2*X*X 
jy=100 . * (y-WMIN) /AMP 
Ny(MX)=i5o-jy 
Ny(MX)=-Ny(MX) 

DO  50  J=2,N1 

x=o. 

A1=WP(J) 

A2=3 . * (W(J+1) -W(J) -WP(J) ) -WP(J+1)+WP(J) 
A3=WP(J+1)-WP(J)-2.*(W(J+1)-W(J)-WP(J) ) 

; CUBIC  FOR  REMAINING,  MATCH 
BOTH  VALUES  AND  SLOPES 


DO  40  1=1, NX 

X=X+DX 

MX=MX+1 

y=W(J)+Al*X+A2*X*X+A3*(X**3) 
jy=100 . * (y-WMIN) /AMP 
Ny(MX)=i50-jy 
Ny(MX)=-Ny(MX) 

CONTINUE 

IF(NFIL.GT.l)  ,v 

OPEN ( 4 , FILE= ' GRAPH . FIL' , FORM= ' FORMATTED ' , STATUS=  NEW  ) 

WRITE (4,100) TITLEl , TITLE2 
WRITE (4,150) NRUN , VALUE 
WRITE(4,200)NORG,MX, (Ny (I) , 1=1 ,MX) 

IF(NFIL.EQ. 1.0R.NFIL.EQ.2)RETURN 
CLOSE ( 4 , STATUS= ' KEEP ' ) 

FORMAT (A1 2) 

FORMAT(I6,E15.6) 

FORMAT (14) 

RETURN 

END 

PROGRAM  PLTGRF 
DIMENSION  W(30) 

CHARACTER  IOUT*72 ,TITLE1*12 ,TITLE2*12 
NFIL=3 

IOUT='  Please  key  in  a  plot  ID  number' 

CALL  MESSAGED tOUT) 

READ (*,10) NRUN 

IOUT='  How  many  data  points?' 

CALL  MESSAGE (lOUT) 

READ(*,10)N 

FORMAT(I8) 

IOUT='  Key-in  each  number,  followed  by  carriage  return 


»•* 


l^'r 


m 


o  o  n  o  o  o  o 


VTnivvwwuiruiftfKvirvTfVTriJrvjrujfVjrujniT 


20 

30 

40 


10 

C 


C 

10 

C 


NADC  88014^ 


CALL  MESSAGE (lOUT) 

READ(*,20) (W(I) ,I*1,N) 

FORMAT (G17. 8) 

IOUT='  Here  are  the  data  points  —  • 

CALL  MESSAGE (lOUT) 

WRITE(*,30) (W(I) ,I=1,N) 

FORMAT (5G15. 6) 

lOUT*'  Please  key-in  two  titles,  each  upto  12  characters' 

CALL  MESSAGE (lOUT) 

READ (*, 4 0 ) TITLE 1 , TITLE2 
FORMAT (A12) 

IOUT='  Please  key-in  a  value  to  go  with  the  last  title' 

CALL  MESSAGE (lOUT) 

READ (*,20) VALUE 

IOUT='  That  is  it.  You  wanna  plot  another?  ' 

CALL  MESSAGE (lOUT) 

IOUT-'  If  yes,  hit  any  NUMBER  key.  If  not,  hit  RETURN' 

CALL  MESSAGE (lOUT) 

READ (*,10) KEY 

IF (NFIL . EQ . 1 . AND . KEY . EQ . 0 ) NFIL*0 
IF (NFIL . EQ . 2 . AND . KEY . NE . 0) NFIL=1 
IF(NFIL.EQ.3.AND.KEY.NE.0)NFIL«2 
IOUT='  Please  wait.  I,in  working  on  the  plot.' 

CALL  MESSAGE (lOUT) 

CALL  GRAPH (W, N , TITLE 1 , TITLE2 , NRUN , VALUE , NFIL) 

IF(KEY.NE.O)GOTO  5 

STOP 

END 

SUBROUTINE  MESSAGE (lOUT) 

CHARACTER  IOUT*72 
WRITE (*, 10 )IOUT 
FORMAT (/A7 2) 

RETURN 

END 

SUBROUTINE  GRAPH ( W , N , TITLEl , TITLE2 , NRUN , VALUE , NFIL) 

SMOOTH  CURVE  FOR  N  EQUALLY  SPACED  POINTS 
NFIL  =  3  FOR  ONE  PLOT  ONLY 

NFIL  =  2  FOR  FIRST,  WITH  MORE  PLOTS  TO  COME 
NFIL  =  1  FOR  MORE  PLOTS  TO  COME 
NFIL  =  0  FOR  LAST 

DIMENSION  W(30),WP(30),NY(640) 

CHARACTER  TITLE1*12 , TITLE2*12 

N1=N-1 
DO  10  1*2, N1 

WP(I)*0.5*(W(I+1)-W(I-1) ) 

WP(N)=0. 

; COERCE  SYMMETRY  ABOUT  LAST  POIN 

WMAX=-l.E+25 
WMIN=l.E+25 
DO  20  1*1, N 
WMAX=AMAX1 (WMAX , W ( I ) ) 
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20 


30 


40 

50 


100 

150 

200 


WMIN«AMIN1 (WMIN , W ( I ) ) 

AMP-WMAX-WMIN 

IF ( AMP. EQ.O.) AMP-1. 

ORG-100. *MMIN/AMP 
NORG-150 .  -t-ORG 
NX-639/N1 
MX-0 

A1-2.*(W(2)-W(1)  )-WP(2) 

A2=WP(2)-W(2)+W(1) 

.•PARABOLIC  FOR  FIRST  INTERVAL, 
MATCH  SLOPE  AT  2ND  POINT 

DX=1./ FLOAT (NX) 

X— DX 

DO  30  I«1,NX 

X-X+DX 

MX-MX+1 

Y-W ( 1 ) +A1*X+A2  *X*X 
JY-100 . * ( Y-WMIN) /AMP 
NY(MX)=150-JY 
NY  (MX)— NY  (MX) 

DO  50  J=2,N1 
X=0. 

A1=WP(J) 

A2=3 . * (W(J+1) -W(J) -WP(J) ) -WP(J+1)+WP(J) 
A3=WP(J+1)-WP(J)-2.*(W(J+1)-W(J)-WP(J) ) 

; CUBIC  FOR  REMAINING,  MATCH 
BOTH  VALUES  AND  SLOPES 

DO  40  1=1, NX 

X=X+DX 

MX=MX+1 

Y=W(J)+A1*X+A2*X*X+A3*(X**3) 

JY-100 . * (Y-WMIN) /AMP 

NY(MX)=150-JY 

NY(MX)=-NY(MX) 

CONTINUE 

IF(NFIL.GT.l) 

OPEN ( 4 , FILE- 'GRAPH . FIL' , FORM- ' FORMATTED ' , STATUS- ' NEW ' ) 
WRITE(4,100)TITLE1,TITLE2 
WRITE (4,150) NRUN , VALUE 
WRITE(4,200)NORG,MX, (NY(I) ,1=1, MX) 
IF(NFIL.EQ.1.0R.NFIL.EQ.2)RETURN 
CLOSE ( 4 , STATUS- ' KEEP ' ) 

FORMAT (A12) 

FORMAT(I6,E15.6) 

FORMAT (14) 

RETURN 

END 
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DISTEIBOTIOH  LIST 

HON-COVEMMtaarr  activities  (continued) 

NO.  OF 
COPIES 

NORTHROP  AIRCRAFT  CORP.,  One  Northrop  Avenue,  Hawthorne,  CA  90250 

(Attn:  Dr.  M.  Ratwanl,  B.  Butler  and  R.  Whitehead).  .  .  3 

PENNSYLVANIA  STATE  UNIVERSITY,  227  Hammond  Building,  Department  of 
Engineering  Science  and  Mechanics,  University  Park,  PA  16802 
(Attn:  H.  Thomas  Hahn,  Ph.D.).  .....  1 

PURDUE  University,  school  of  Aeronautics  and  Astronaut '.cs. 

West  Lafayette,  IN  47907 

(Attn:  Dr.  C.  T.  Sun).  ....  ...1 

ROCKWELL  INTERNATIONAL,  Columbus.  OH  43216 

(Attn:  M.  Schvelger).  .......  1 

ROCKWELL  INTERNATIONAL,  Los  Angeles,  CA  90009 

(Attn:  Dr.  Lackman).  .......  1 

(Attn:  W.  O'Brien).  .......  1 

ROCKWELL  INTERNATIONAL,  Tulsa,  OK  74151 

(Attn:  F.  Kaufman).  .......  1 

SIKORSKY  AIRCRAFT,  Stratford,  CT  06622 

(Attn:  S.  Garbo).  .......1 

J.  P.  STEVENS  &  CO.,  INC.,  New  York,  NY  10036 

(Attn:  H.  I.  Shulock).  .......  1 

TUSKEGEE  UNIVERSITY,  School  of  Engineering  and  Architecture, 

Tuskegee,  AL  36088 

(Attn:  Vascar  G.  Harris,  Dean).  .....  1 

UNIVERSITY  OF  DAYTON  RESEARCH  INSTITUTE,  300  College  Park  Avenue, 

Dayton,  OH  45469 

(Attn:  Dr.  J.  Gallagher).  ......  1 

UNIVERSITY  OF  DELAWARE,  Mechanics  &  Aerospace  Eng.  Dept., 

Evans  Hall,  Newark,  DE  19711 

(Attn:  Dr.  R.  B.  Pipes,  Dr.  J.  R.  Vinson  and  Dr.  D.Wilkins  .  3 

UNIVERSITY  OF  OKLAHOMA,  Norman,  OK  73019 

(Attn:  Dr.  C.  W.  Bert,  School  of  AMNE).  ....  1 

UNIVERSITY  OF  WYOMING,  Laramie,  WY  82071 

(Attn:  Dr.  D.  Adams).  .......  1 

VILLANOVA  UNIVERSITY,  Villanova,  PA  19085 

(Attn:  Dr.  P.  V.  McLaughlin).  .....  1 

VIRGINIA  POLYTECHNIC  INSTITUTE,  Blacksburg,  VA  24061 

(Attn:  Dr.  K.  Relfsnlder).  ......  1 
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PISTRIBOnON  LIST 

NON-GOVERNMENT  ACTIVITIKS  (continued) 

NO.  OF 
COPIES 

GENERAL  ELECTRIC  CO.,  Philadelphia,  PA  19101 

(Attn:  A.  Garber,  C.  Zweben).  .....  2 

GRUMMAN  CORPORATION,  South  Oyster  Bay  Rd.,  Bethpage,  NY  11714 

(Attn:  R.  Radcoclc)  ........1 

(Attn:  S.  Pasting.  .......1 

HITCO,  1600  West  13Sth  Street,  Gardena,  CA  90249 

(Attn:  N.  Myers).  .......1 

ITT  RES^CH  INSTITUTE,  Chicago,  IL  60616 

(Attn:  K.  Hofar).  .......1 

KAMAN  AIRCRAFT  CORP.,  Bloomfield,  CT  06002 

(Attn:  Technical  Library).  ......  1 

LEHIGH  UNIVERSITY,  Bethlehem,  PA  18015 

(Attn:  Dr.  G.  C.  Slh).  .......1 

LEONARD  ASSOCIATES,  INC.,  6  East  Avenue,  Mt.  Carmel,  PA  17851 

(Attn:  Mr.  L.  Marchlnskl).  ......  1 

LOCKHEED-CALIFORNIA  CO.,  Burbank,  CA  91520 

(Attn:  E.  K.  Walker )..  .  .  .  .  .  .  1 

(Attn:  A.  «James).  .......1 

LOCKHEED-MISSILES  &  SPACE  CO.,  1111  Lockheed  Way,  Sunnyvale,  CA  94086 

(Attn:  .  A.  Bailie).  .......1 

LOCKHEED-CALIFORNIA  CO.,  Rye  Canyon  Research  Laboratory, 

Burbank,  CA  91520 

(Attn:  D.  E.  Pettit).  .......1 

LOCKHEED-GEORGIA  CO.,  Marietta,  GA  30063 

(Attn:  Technical  Information  Dept.,  72-34,  Zone  26).  .  .  1 

LTV  AEROSPACE  &  DEFENSE  CO. ,  Vought  Missile  &  Advanced  Program 
Division,  P.O.  Box  225907,  Dallas,  TX  75265-0003 

(Attn:  R.  Knight).  .......1 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY,  Technology  Laboratory  for 

Advanced  Composite,  77  Massachuetts  Avenue,  Cambridge,  MA  02139 
(Attn:  Dr.  P.  A.  Lagace).  ......  1 

MATERIALS  SCIENCES  CORP.,  Spring  House,  PA  19477 

(Attn:  Dr.  B.  W.  Rosen).  ......  1 

McDONNELL-DOUGLAS  CORP.,  St.  Louis,  MO  63166 

(Attn:  R.  Plnckert.  .  .  .  .  .  .  .  .  1 

McDONNELL-DOUGLAS  CORP.  Long  Beach,  CA  90846 

(Attn:  C.  Y.  Kam,  Dept.  C1-E84,  MC  212-10).  ...  1 

McDONNELL-DOUGLAS  HELICOPTER  CO.,  Bldg.  530,  5000  E.  McDowell  Rd 
Mesa,  AZ  85258 

(Attn:  J.  K.  Sen,  M.S.  338).  ......  1 

McDONNELL-DOUGLAS  HELICOPTER  CO.,  5000E.  McDowell,  M/S  B337 
Mesa,  AZ  85205 
(Attn:  Steve  Guymon). 
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DISTRIBDTIOH  LIST 
MOH-GOVKRMMENT  ACTIVITIES 

NO.  OF 
COPIES 

ANAM^-T  LABORATORIES,  3400  Investment  Blvd.,  Hayward,  CA  94545-3811 

(Attn:  Dr.  R.  Arnold).  .......  1 

ALCOA  DEFENSE  SYSTEMS  CORP.,  16761  Via  delCampo  Court, 

San  Diego,  CA  92127 

(Attn:  D.  Myers).  .......  1 

AVCO,  Specialty  Materials  Dlv. ,  2  Industrial  Avenue,  Lovell,  MA  01851 

(Attn:  Mr.  W.  F.  Grant).  ......  1 

BATTELLE  COLUMBUS  LABORATORIES,  Metals  and  Ceramics  Information  Center 

505  King  Avenue,  Columbus,  OU  43201.  ....  1 

BEECH  AIRCRAFT  CORP.,  4130  Linden  Avenue,  Dayton,  OH  45432 

(Attn:  M.  B.  Goetz).  .......  1 

BELL  AEROSPACE  COMPANY,  Buffalo,  NY  14240 

(Attn:  F.  M.  Anthony,  Zone  1-85).  .....  1 

BELL  HELICOPTER  CO.,  Fort  Worth,  TX  76101 

(Attn:  M.  K.  Stevenson).  ......  1 

BENDIX  PRODUCTS,  Aerospace  Division,  South  Bend,  IN  46619 

(Attn:  R.  V.  Cervelll).  ......  1 

BOEING  CO.,  P.  0.  Box  3707,  Seattle,  WA  98124 

(Attn:  J.  Qulnllven,  and  Dr.  R.  June).  ....  2 

BOEING  HELICOPTER  COMPANY,  P.O.  Box  16858,  Philadelphia,  PA  19143 

(Attn:  R.  L.  Pinckney).  ......  1 

(Attn:  D.  Hart).  ........  1 

(Attn:  C .  Albrecht )..  .  .  .  .  .  .  1 

BOEING  CO.,  Wichita,  KS  67277-7730.  .  1 

(Attn:  J.  Avery ).  .  .  .  .  .  .  .  1 

(Attn:  R.  Waner). .  .......1 

CABOT  CORPORATION,  Billerica  Research  Center,  Billerica,  MA  01821.  1 

DEPARTMENT  OF  TRANSPORTATION,  Kendall  Square,  Cambridge,  MA  02142.  1 

(Attn:  Dr.  Ping  Tong,  DTS  76,  TSC).  ....  1 

DREXEL  UNIVERSITY,  Philadelphia,  PA  19104 

(Attn:  Dr.  P.  C.  Chou).  ......  1 

(Attn:  Dr.  A.  S.  D.  Wang).  ......  1 

E.  I.  DuPONT  COMPANY,  Textile  Fibers  Department,  Chestnut  Run  Location 
CR701,  Wilmington,  DE  19898 

(Attn:  V.  L.  Bertarelli).  ......  1 

GEORGIA  INSTITUTE  CF  TECHNOLOGY,  Atlanta,  GA  30332 

(Attn:  (L.  Fehfield)).  .......  1 

GENERAL  DYNAMICS /CO WAIR,  San  Diego,  CA  92138 

(Attn:  Dr.  R  Dunbar).  .......  1 

GENERAL  DYNAMICS,  Fi.rt  Worth  Division,  PO  Box  748,  Fort  Worth,  TX  76101 
(Attn:  J.  A.  F ant )..  .  .  .  .  .  .  I 

(Attn:  Composite  Structures  Eng.  Dept.).  ....  1 
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NAVSHIPRANDCEN,  Annapolis,  MD  21403 

(Attn:  H.  Edlesteln,  Code  2870).  .... 

NRL,  Washington,  D.C.  20375 

(Attn:  Dr.  I.  Wolock,  Code  6122;  Dr.  C.  I.  Chang, 

and  Dr.  R.  Badallance).  .... 

NSWC,  WHITE  OAK  LABORATORY,  Silver  Spring,  MD  20910 

(Attn:  Dr.  J.  Goff,  Materials  Evaluation  Branch,  Code  R-34 
(Attn:  Dr.  JT.  M.  Augl ) .  .  .  .  .  . 

ONR,  800  N.  Quincy  St.,  Arlington,  VA  22217 

(Attn:  A.  Kushner  Code  1132SM;  Y.  Rajapakse,  Code  1132SM) 
(Attn:  R.  Jones,  Code  1216).  ..... 

ONT,  800  N.  Quincy  Street,  Arlington,  VA  22217 

(Attn:  Capt.  K.  Cox,  (0NT-21D).  .... 

PLASTEC,  Plcatinny  Arsenal,  Dover,  NJ  07801 

(Attn:  H.  Pebly).  ...... 

(Attn:  Librarian,  Code  DRDAR-SCM-0,  Bldg.  351-N). 

ARMY  MATERIALS  TECHNOLOGY  LABORATORY 
Watertown,  MA  02171 

(Attn:  D.  Oplinger ,  SLCMT-MS).  .... 

U.  S.  ARMY  APPLIED  TECHNOLOGY  LABORATORY,  USARTL,  (AVRADCOM), 

Ft.  Eustis,  VA  23604-5418 

(Attn:  J.  Waller;  G.  McAllister.  .... 

U.  S.  ARMY  AIR  MOBILITY  R&D  LABORATORY,  Ft.  Eustis,  VA  23604 
(Attn:  H.  Reddick).  ...... 

U.  S.  ARMY  R&T  LABORATORY  (AVRADCOM),  Ames  Research  Center, 

Moffet  Field,  CA  94035 

(Attn:  F.  Immen,  DAVDL-AS-MS  207-5). 

U.  S.  NAVAL  ACADEMY,  Annapolis,  MD  21402 

(Attn:  Mechanical  Engineering  Department). 

DAVID  TAYLOR  NAVAL  SHIP  RESEARCH  &  DEVELOPMENT  CENTER, 

Annapolis,  MD  21402 

(Attn:  E.  T.  Camponeschi,  Code  2844;  R.  Crane,  Code  2844). 
DAVID  TAYLOR  NAVAL  SHIP  R&D  CENTER 
Bethesda,  MD  20084 

(Attn:  A.  Macander,  Code  1720).  .... 


NAVAIRDEVCEN,  Warminster,  PA  18974 
(Attn:  Code  8131). 

(Attn:  Code  09L2). 
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AFWAL,  WPAFB,  OH  45433-6553 

(Attn:  FIBEC,  Dr.  G.  SendeckyJ).  ..... 

(Attn:  FIBC/L.  Kelly) . 

(Attn:  FIBCA/C.  Ramsey).  ...... 

(Attn:  FIBG/H.  F.  Wolff) . 

AFWAL,  WPAFB,  OH  45433-6533 

(Attn:  MLBM/Dr.  J.  Whitney,  M.  Knight).  .... 

(Attn:  MLB/F.  Cherry).  ....... 

(Attn:  MBC/Relnhart ) .  ....... 

(Attn:  MLSE/S.  Fecheck).  ...... 

DEPARTMENT  OF  THE  AIR  FORCE,  Bldg.  410,  Bolling  Air  Force  Base, 
Washington,  D.C.  20332 

(Attn:  Dr.  M.  Salklnd,  Dr.  Amos).  ..... 
DEFENSE  TECHNICAL  INFORMATION  CENTER  (DTIC),  Bldg. #5,  Cameron  Station 
Alexandria,  VA  22314 

(Attn:  Administrator).  ....... 

FAA,  Washington,  D.C.  20591 

(Attn:  J.  R.  Soderqulst,  AW-103).  ..... 

FAA,  Technical  Center,  Atlantic  City,  NJ  08405 

(Attn:  L.  Neri ,  Code  ACT-330;  C.  Caiafa,  Code  ACT-033). 

NASA  HEADQUARTERS,  Washington,  D.  C.  20546 

(Attn:  Airframes  Branch,  FS-120).  ..... 

(Attn:  OAST/RM,  Dr.  D.  Mulvllle).  ..... 

NASA,  George  C.  Marshall  Space  Flight  Center,  Huntsville,  AL  35812 
(Attn:  R.  Schwinghamer ,  S&E-ASTN-M).  .... 

NASA,  Langley  Research  Center,  Hampton,  VA  23365 

(Attn:  Dr.  J.  R.  Davidson,  MS  188E;  Dr.  J.  Starnes,  MS-190; 

Dr.  M.  Mikulus,  H.  Bohan,  and  Dr.  C.  P.  Blakenshlp 
MS  189M).  ........ 

NASA,  Lewis  Research  Center,  Cleveland,  OH  44135 

(Attn:  Dr.  C.  Chamis,  MS  49-6;  M.  Hershberg,  MS  49-6). 
NAVAIRSYSCOM,  Washington,  D.C.  20361 

(Attn:  AIR-00D4).  ....... 

(Attn:  AIR-530).  ........ 


(Attn: 

(Attn:  AIR-5302D). 

(Attn:  AIR-5302). 

(Attn:  AIR-5302F). 

(Attn:  AIR-53032D). 

(Attn:  AIR-931B). 

NAVPGSCHL,  Monterey,  CA  95940 

(Attn:  Prof.  R.  Ball,  Prof.  M.  H. 
NAVSEASYSCOM,  Washington,  D.C.  20360 
(Attn:  C.  Zannls,  SEA-05R25)). 
NSEC,  Arlington,  VA  20360 
(Attn:  NSEC-6101E). 


Bank,  Prof.  K.  Challenger). 
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